Proofs & evidence

{
  "schema": "rheon-bounded-integration-qualification-v1",
  "status": "checked",
  "recorded_utc": "2026-10-04T11:30:12.015392+00:00",
  "qualified_source_commit": "c6809318d8973dc9a2e3af130eac1515c68cddb3",
  "qualified_source_tree": "90572bf83a66fa5c93a7f9fb763e4aa071ce8a74",
  "integration_base": "2bc2a84fd4c562062ad592426794e8382f098a77",
  "bounded_source_commit": "0ba3135754c014ed1e08e945601c401823835a5a",
  "bounded_evidence_commit": "6463d71476dbab4aacd06094a96d8a5f567d2f2f",
  "bounded_mapped_commits": [
    "bdb8320a8b5d160a41d392db17a94432ddc49e2d",
    "eaba5ab960fa1e895b5b225fb2f976c5c2941574"
  ],
  "public_theorems": 42,
  "previous_public_theorems": 33,
  "new_public_theorems": 9,
  "new_definitions": 6,
  "audited_declarations": 60,
  "focused_audited_declarations": 15,
  "allowed_transitive_axioms": [
    "propext",
    "Classical.choice",
    "Quot.sound"
  ],
  "lean_version": "4.19.0",
  "lean_release_sha256": "6fe3ce97a58f44e2b3567d455b994eacec5bfe9ae7774f2a573444480ba813fe",
  "dependencies": [
    {
      "name": "mathlib",
      "revision": "c44e0c8ee63ca166450922a373c7409c5d26b00b",
      "tracked_sources_match": true
    },
    {
      "name": "plausible",
      "revision": "77e08eddc486491d7b9e470926b3dbe50319451a",
      "tracked_sources_match": true
    },
    {
      "name": "LeanSearchClient",
      "revision": "25078369972d295301f5a1e53c3e5850cf6d9d4c",
      "tracked_sources_match": true
    },
    {
      "name": "importGraph",
      "revision": "e6a9f0f5ee3ccf7443a0070f92b62f8db12ae82b",
      "tracked_sources_match": true
    },
    {
      "name": "proofwidgets",
      "revision": "c4919189477c3221e6a204008998b0d724f49904",
      "tracked_sources_match": true
    },
    {
      "name": "aesop",
      "revision": "5d50b08dedd7d69b3d9b3176e0d58a23af228884",
      "tracked_sources_match": true
    },
    {
      "name": "Qq",
      "revision": "fa4f7f15d97591a9cf3aa7724ba371c7fc6dda02",
      "tracked_sources_match": true
    },
    {
      "name": "batteries",
      "revision": "f5d04a9c4973d401c8c92500711518f7c656f034",
      "tracked_sources_match": true
    },
    {
      "name": "Cli",
      "revision": "02dbd02bc00ec4916e99b04b2245b30200e200d0",
      "tracked_sources_match": true
    }
  ],
  "method": "Fresh native root module build and shared/focused audits with existing pinned dependency build cache; no dependency/pin changes or network workaround. Source gate and four rejection fixtures pass normally and under optimized Python; negative custom/admitted/skipped-name probes pass in both modes.",
  "source_inventory": {
    "AxiomAudit.lean": "ff82b9cc304ed6897fb7e7c252fdae9988b86d3c060f66924f3ec7f136054bbb",
    "Rheon/BoundedPhysics.lean": "20dd53b1b137a01f8b04de47ffbe8474561696940e17ab1ed4d3bfde388b03af",
    "Rheon/Discrete.lean": "70a2ca5d2c5cbefacbbb8742954375fdfe8cc25c3c5b8afd05738941cfc5ee0a",
    "Rheon/Indexing.lean": "862c33164e75daccd43603c9b5a61473fd825885d9106b766a35204e578f60bf",
    "Rheon/Physics.lean": "26cbb4f75076c7e70a90da126f96622b75044befce26a4c17a50066c76e96050",
    "Rheon/Transport.lean": "f967e593db1df0ec945ea6830e67b991e56da168e9f6e3f2d42c4c5866c5cf06",
    "Rheon.lean": "94211dc1dbf7cadb5a39f441aecd690c4c619550e15e3fedf5b6a63ff3ce86a3"
  },
  "bounded_evidence_file_sha256": {
    "proofs/Rheon/BoundedPhysics.lean": "20dd53b1b137a01f8b04de47ffbe8474561696940e17ab1ed4d3bfde388b03af",
    "evidence/bounded-physics-contracts/AuditBoundedPhysics.lean": "ad9b42dd49701214c1d90d8836eccbc5ea4ad1ea06ffdeee52748d7e190ce647",
    "evidence/bounded-physics-contracts/README.md": "28826e941006c4169da29e09e8e072fbb92b7238952a7fd42250aefa5268ea3a",
    "evidence/bounded-physics-contracts/audit-initial-compiler-body-failure.log": "360b9115329f1ff722f2e1139a6275b445befe045f86a849df9d8fb2230dbbc8",
    "evidence/bounded-physics-contracts/audit.log": "f3f431c2f6aacb156cb3045bc852b7c035e9ad049984ff42c7671ef3bd97bf61",
    "evidence/bounded-physics-contracts/build-initial-warning.log": "230f69f153e1bb930275577ff1bbc75144aacab7eec449330f62ad6f690c5b45",
    "evidence/bounded-physics-contracts/build.log": "0ea66f7f3cc115a8f2c378e6ca34cdcba10b3d6e26faf66844c27c5203b2a57d",
    "evidence/bounded-physics-contracts/prerequisite-build.log": "dcbc06c6703ef8a53d4ed5f8ffc38b92810368868e0869c733f41811197a5470",
    "evidence/bounded-physics-contracts/witness-optimized.json": "b1cfd020eec82b1e67fca6dba3a4dc40dc67b76722c59038a78ec299e8270846",
    "evidence/bounded-physics-contracts/witness.json": "b1cfd020eec82b1e67fca6dba3a4dc40dc67b76722c59038a78ec299e8270846",
    "evidence/bounded-physics-contracts/witness.py": "b7cd7fb53b855a53a9f71f1b734ae3a09fbf523ea358f3cef6860f4fd7534df6"
  },
  "logs": {
    "witness-optimized.json": "b1cfd020eec82b1e67fca6dba3a4dc40dc67b76722c59038a78ec299e8270846",
    "focused-audit.log": "f3f431c2f6aacb156cb3045bc852b7c035e9ad049984ff42c7671ef3bd97bf61",
    "lean-build.log": "46a30e4abc68b8894fd95b721dae9fa364c65a93e267c140fa981a41525f35af",
    "axiom-audit.log": "4e8c08d56c73e1bec1f359636a53454a6cdca66a9288415511bb3ce2c1b59371",
    "witness.json": "b1cfd020eec82b1e67fca6dba3a4dc40dc67b76722c59038a78ec299e8270846",
    "negative-audit.log": "6c31afa942a90b8f697f122ef8f310c31ea794a3aacffb7dcf0cacd8724a7c52"
  },
  "rational_witnesses": "normal and optimized outputs byte-identical; crossing/clipping, endpoint-sign omissions, positive/negative viscosity and omitted-update examples",
  "limits": "Strict crossing of an infinite stationary plane; arbitrary fixed finite strain with exact unforced homogeneous coordinate update and nonnegative weights/masses/time. No facet containment, earliest arbitrary-mesh query, assembled symmetric-gradient/traction, solve existence/uniqueness, forcing, approximate solve, PDE or IEEE refinement."
}

The pinned reconstruction passed its normal lake build, allowed-axiom audit and rejection/source-policy probes. Those counts and receipts apply to its recorded source. Earlier PR24 57/77 qualification and recorded native lab identities remain historical.

Associated reconstruction Lean CI · Source identities, receipts and historical bindings

Download Lean source

import Rheon.BoundedPhysics
import Rheon.ObstacleOperators
import Mathlib.Algebra.Order.BigOperators.Ring.Finset
import Mathlib.Data.Finset.Fold

/-! Exact finite-row contracts for the reconstructed aligned strain operator.
The coefficient bound is derived, including zero rows, from the per-face
absolute-row bound. Assembly, represented geometry, IEEE arithmetic and
approximate projection are external obligations. -/
namespace Rheon.AlignedStrain
noncomputable section
open scoped BigOperators
open Rheon.ObstacleOperators

def gather {n r : ℕ} (E : Fin r → Fin n → ℝ) (u : Fin n → ℝ) (q : Fin r) : ℝ :=
  ∑ f, E q f * u f

def strainOperator {n r : ℕ} (E : Fin r → Fin n → ℝ) (w : Fin r → ℝ)
    (u : Fin n → ℝ) (f : Fin n) : ℝ :=
  ∑ q, E q f * w q * gather E u q

def strainLoss {n r : ℕ} (E : Fin r → Fin n → ℝ) (w : Fin r → ℝ)
    (u : Fin n → ℝ) : ℝ := ∑ q, w q * (gather E u q)^2

def rowAbs {n r : ℕ} (E : Fin r → Fin n → ℝ) (q : Fin r) : ℝ :=
  ∑ f, |E q f|

def faceBound {n r : ℕ} (E : Fin r → Fin n → ℝ) (w : Fin r → ℝ)
    (mass : Fin n → ℝ) (f : Fin n) : ℝ :=
  (∑ q, w q * |E q f| * rowAbs E q) / mass f

/-- Maximum of face bounds, with zero for the empty active space. -/
def rowBound {n r : ℕ} (E : Fin r → Fin n → ℝ) (w : Fin r → ℝ)
    (mass : Fin n → ℝ) : ℝ := Finset.univ.fold max 0 (faceBound E w mass)

theorem row_bound_nonnegative {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (mass : Fin n → ℝ) : 0 ≤ rowBound E w mass := by
  exact (Finset.le_fold_max 0).mpr (Or.inl le_rfl)

theorem face_le_row_bound {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (mass : Fin n → ℝ) (f : Fin n) :
    faceBound E w mass f ≤ rowBound E w mass := by
  exact (Finset.le_fold_max (faceBound E w mass f)).mpr
    (Or.inr ⟨f, Finset.mem_univ f, le_rfl⟩)

def forceNorm {n r : ℕ} (E : Fin r → Fin n → ℝ) (w : Fin r → ℝ)
    (mass u : Fin n → ℝ) : ℝ :=
  ∑ f, (strainOperator E w u f)^2 / mass f

def euler {n r : ℕ} (E : Fin r → Fin n → ℝ) (w : Fin r → ℝ)
    (mass u : Fin n → ℝ) (dt mu : ℝ) (f : Fin n) : ℝ :=
  u f - dt * mu * strainOperator E w u f / mass f

theorem transpose_work {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (u : Fin n → ℝ) :
    (∑ f, strainOperator E w u f * u f) = strainLoss E w u := by
  exact Rheon.BoundedPhysics.viscous_work E w u

theorem loss_nonnegative {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (u : Fin n → ℝ) (hw : ∀ q, 0 ≤ w q) :
    0 ≤ strainLoss E w u := by
  exact Finset.sum_nonneg (fun q _ => mul_nonneg (hw q) (sq_nonneg _))

theorem transpose_pairing {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (u v : Fin n → ℝ) :
    (∑ f, strainOperator E w u f * v f) =
      ∑ q, w q * gather E u q * gather E v q := by
  unfold strainOperator
  simp only [Finset.sum_mul]
  rw [Finset.sum_comm]
  apply Finset.sum_congr rfl
  intro q _
  calc
    (∑ f, E q f * w q * gather E u q * v f) =
        ∑ f, (w q * gather E u q) * (E q f * v f) := by
      apply Finset.sum_congr rfl
      intro f _
      ring
    _ = w q * gather E u q * gather E v q := by
      rw [← Finset.mul_sum]
      rfl

theorem operator_symmetric {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (u v : Fin n → ℝ) :
    (∑ f, strainOperator E w u f * v f) =
      ∑ f, strainOperator E w v f * u f := by
  rw [transpose_pairing, transpose_pairing]
  apply Finset.sum_congr rfl
  intro q _
  ring

theorem row_abs_nonnegative {n r : ℕ} (E : Fin r → Fin n → ℝ) (q : Fin r) :
    0 ≤ rowAbs E q := Finset.sum_nonneg (fun _f _ => abs_nonneg _)

theorem coefficient_le_row_abs {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (q : Fin r) (f : Fin n) : |E q f| ≤ rowAbs E q := by
  exact Finset.single_le_sum (fun g _ => abs_nonneg (E q g)) (Finset.mem_univ f)

theorem zero_row_coefficient {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (q : Fin r) (h : rowAbs E q = 0) (f : Fin n) : E q f = 0 := by
  have hb := coefficient_le_row_abs E q f
  rw [h] at hb
  exact abs_eq_zero.mp (le_antisymm hb (abs_nonneg _))

theorem zero_row_gather {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (u : Fin n → ℝ) (q : Fin r) (h : rowAbs E q = 0) : gather E u q = 0 := by
  unfold gather
  simp_rw [zero_row_coefficient E q h, zero_mul]
  simp

/-- Weighted Cauchy at one face; the zero-row branch needs no division. -/
theorem face_cauchy {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (u : Fin n → ℝ) (hw : ∀ q, 0 ≤ w q) (f : Fin n) :
    (strainOperator E w u f)^2 ≤
      (∑ q, w q * |E q f| * rowAbs E q) *
        ∑ q, w q * |E q f| * (gather E u q)^2 / rowAbs E q := by
  apply Finset.sum_sq_le_sum_mul_sum_of_sq_eq_mul Finset.univ
    (r := fun q => E q f * w q * gather E u q)
    (f := fun q => w q * |E q f| * rowAbs E q)
    (g := fun q => w q * |E q f| * (gather E u q)^2 / rowAbs E q)
  · intro q _
    exact mul_nonneg (mul_nonneg (hw q) (abs_nonneg _)) (row_abs_nonnegative E q)
  · intro q _
    exact div_nonneg (mul_nonneg (mul_nonneg (hw q) (abs_nonneg _)) (sq_nonneg _))
      (row_abs_nonnegative E q)
  · intro q _
    by_cases h : rowAbs E q = 0
    · simp [h, zero_row_coefficient E q h f]
    · field_simp [h]
      ring_nf
      simp only [sq_abs]
      ring

/-- The per-face coefficient bound proves K M^-1 K <= B K as quadratic forms. -/
theorem coefficient_force_bound {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (mass u : Fin n → ℝ) (B : ℝ)
    (hw : ∀ q, 0 ≤ w q) (hm : ∀ f, 0 < mass f)
    (_hB : 0 ≤ B) (bound : ∀ f, faceBound E w mass f ≤ B) :
    forceNorm E w mass u ≤ B * strainLoss E w u := by
  have hf : ∀ f, (strainOperator E w u f)^2 / mass f ≤
      B * ∑ q, w q * |E q f| * (gather E u q)^2 / rowAbs E q := by
    intro f
    have hn : 0 ≤ ∑ q, w q * |E q f| * (gather E u q)^2 / rowAbs E q := by
      apply Finset.sum_nonneg
      intro q _
      exact div_nonneg (mul_nonneg (mul_nonneg (hw q) (abs_nonneg _)) (sq_nonneg _))
        (row_abs_nonnegative E q)
    have hb : (∑ q, w q * |E q f| * rowAbs E q) ≤ B * mass f :=
      (div_le_iff₀ (hm f)).mp (bound f)
    apply (div_le_iff₀ (hm f)).mpr
    calc
      (strainOperator E w u f)^2 ≤ _ := face_cauchy E w u hw f
      _ ≤ (B * mass f) * _ := mul_le_mul_of_nonneg_right hb hn
      _ = _ := by ring
  have hs : (∑ f, ∑ q, w q * |E q f| * (gather E u q)^2 / rowAbs E q) =
      strainLoss E w u := by
    rw [Finset.sum_comm]
    unfold strainLoss
    apply Finset.sum_congr rfl
    intro q _
    by_cases h : rowAbs E q = 0
    · simp [h, zero_row_gather E u q h]
    · calc
        (∑ f, w q * |E q f| * (gather E u q)^2 / rowAbs E q) =
            (w q * (gather E u q)^2 / rowAbs E q) * ∑ f, |E q f| := by
          rw [Finset.mul_sum]
          apply Finset.sum_congr rfl
          intro f _
          ring
        _ = w q * (gather E u q)^2 := by
          change (w q * (gather E u q)^2 / rowAbs E q) * rowAbs E q = _
          exact div_mul_cancel₀ _ h
  calc
    forceNorm E w mass u ≤ ∑ f, B *
        ∑ q, w q * |E q f| * (gather E u q)^2 / rowAbs E q :=
      Finset.sum_le_sum (fun f _ => hf f)
    _ = B * strainLoss E w u := by rw [← Finset.mul_sum, hs]

/-- Exact coordinate Euler energy expansion, before any step restriction. -/
theorem euler_energy_identity {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (mass u : Fin n → ℝ) (dt mu : ℝ)
    (hm : ∀ f, mass f ≠ 0) :
    kineticEnergy mass (euler E w mass u dt mu) - kineticEnergy mass u =
      -(dt * mu) * strainLoss E w u + (dt * mu)^2 / 2 * forceNorm E w mass u := by
  have local_identity : ∀ f, mass f / 2 * (u f - dt * mu *
      strainOperator E w u f / mass f)^2 - mass f / 2 * (u f)^2 =
      -(dt * mu) * (strainOperator E w u f * u f) +
        (dt * mu)^2 / 2 * ((strainOperator E w u f)^2 / mass f) := by
    intro f
    field_simp [hm f]
    ring
  unfold kineticEnergy euler
  rw [← Finset.sum_sub_distrib]
  calc
    _ = ∑ f, (-(dt * mu) * (strainOperator E w u f * u f) +
        (dt * mu)^2 / 2 * ((strainOperator E w u f)^2 / mass f)) := by
      apply Finset.sum_congr rfl
      intro f _
      exact local_identity f
    _ = -(dt * mu) * strainLoss E w u +
        (dt * mu)^2 / 2 * forceNorm E w mass u := by
      rw [Finset.sum_add_distrib, ← Finset.mul_sum, ← Finset.mul_sum, transpose_work]
      rfl

/-- Unforced exact-real Euler decrease under the derived coefficient restriction. -/
theorem euler_energy_nonincrease {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (mass u : Fin n → ℝ) (dt mu B : ℝ)
    (hw : ∀ q, 0 ≤ w q) (hm : ∀ f, 0 < mass f)
    (ht : 0 ≤ dt) (hmu : 0 ≤ mu) (hB : 0 ≤ B)
    (bound : ∀ f, faceBound E w mass f ≤ B) (step : dt * mu * B ≤ 2) :
    kineticEnergy mass (euler E w mass u dt mu) ≤ kineticEnergy mass u := by
  have hb := coefficient_force_bound E w mass u B hw hm hB bound
  have hd := loss_nonnegative E w u hw
  have hi := euler_energy_identity E w mass u dt mu (fun f => ne_of_gt (hm f))
  have t0 : 0 ≤ dt * mu := mul_nonneg ht hmu
  have hterm : (dt * mu)^2 / 2 * forceNorm E w mass u ≤
      (dt * mu)^2 / 2 * (B * strainLoss E w u) :=
    mul_le_mul_of_nonneg_left hb (div_nonneg (sq_nonneg _) (by norm_num))
  have hc : dt * mu * B * strainLoss E w u ≤ 2 * strainLoss E w u :=
    mul_le_mul_of_nonneg_right step hd
  have hc' := mul_le_mul_of_nonneg_left hc t0
  nlinarith

/-- The computed finite maximum supplies the bound premises automatically. -/
theorem euler_row_bound_nonincrease {n r : ℕ} (E : Fin r → Fin n → ℝ)
    (w : Fin r → ℝ) (mass u : Fin n → ℝ) (dt mu : ℝ)
    (hw : ∀ q, 0 ≤ w q) (hm : ∀ f, 0 < mass f)
    (ht : 0 ≤ dt) (hmu : 0 ≤ mu) (step : dt * mu * rowBound E w mass ≤ 2) :
    kineticEnergy mass (euler E w mass u dt mu) ≤ kineticEnergy mass u := by
  exact euler_energy_nonincrease E w mass u dt mu (rowBound E w mass)
    hw hm ht hmu (row_bound_nonnegative E w mass) (face_le_row_bound E w mass) step

/-- Composition uses exactly the pressure masses and exact full residual. -/
theorem euler_exact_projection_nonincrease {n r c : ℕ}
    (E : Fin r → Fin n → ℝ) (w : Fin r → ℝ)
    (P : Fin c → Fin n → ℝ) (area distance u : Fin n → ℝ) (p : Fin c → ℝ)
    (rho dt mu B : ℝ) (hw : ∀ q, 0 ≤ w q)
    (hr : 0 < rho) (ha : ∀ f, 0 < area f) (hd : ∀ f, 0 < distance f)
    (ht : 0 < dt) (hmu : 0 ≤ mu) (hB : 0 ≤ B)
    (bound : ∀ f, faceBound E w (faceMass rho area distance) f ≤ B)
    (step : dt * mu * B ≤ 2)
    (solved : ∀ i, pressureResidual P rho dt area distance
      (euler E w (faceMass rho area distance) u dt mu) p i = 0) :
    kineticEnergy (faceMass rho area distance)
      (pressureCorrected P rho dt distance
        (euler E w (faceMass rho area distance) u dt mu) p) ≤
      kineticEnergy (faceMass rho area distance) u := by
  have hm : ∀ f, 0 < faceMass rho area distance f :=
    fun f => mul_pos (mul_pos hr (ha f)) (hd f)
  exact le_trans
    (exact_pressure_energy_nonincrease P rho dt area distance _ p hr ha hd
      (ne_of_gt ht) solved)
    (euler_energy_nonincrease E w _ u dt mu B hw hm (le_of_lt ht) hmu hB bound step)

/-- The reflected unit corner has exactly the claimed coupled quadratic block. -/
theorem unit_corner_block (U V : ℝ) :
    (1 / 4 : ℝ) * (2 * U + 2 * V)^2 + (1 / 4 : ℝ) * (2 * U)^2 +
      (1 / 4 : ℝ) * (2 * V)^2 = 2 * U^2 + 2 * U * V + 2 * V^2 := by ring

/-- Geometric area equality is an explicit extra volume hypothesis. -/
theorem flat_conductance (volume area delta : ℝ) (hd : delta ≠ 0)
    (product : volume = area * delta) : volume / delta^2 = area / delta := by
  rw [product]
  field_simp [hd]
  ring

end
end Rheon.AlignedStrain

Download Lean source

import Rheon.Physics

/-! Exact finite contracts for a planar collision and an unforced viscous step.
They do not certify mesh queries, stencil assembly, material fits or IEEE code.
The wall normal points into the permitted half-space; it need not be unit length.
The strain matrix and viscosity weights are fixed during the implicit step. -/
namespace Rheon.BoundedPhysics
noncomputable section
open scoped BigOperators

def segment (a b : Fin 3 → ℝ) (t : ℝ) (i : Fin 3) : ℝ :=
  (1 - t) * a i + t * b i

def wallValue (normal : Fin 3 → ℝ) (offset : ℝ) (x : Fin 3 → ℝ) : ℝ :=
  (∑ i, normal i * x i) - offset

def wallHit (normal : Fin 3 → ℝ) (offset : ℝ) (a b : Fin 3 → ℝ) : ℝ :=
  wallValue normal offset a / (wallValue normal offset a - wallValue normal offset b)

theorem wall_segment_affine (normal : Fin 3 → ℝ) (offset : ℝ)
    (a b : Fin 3 → ℝ) (t : ℝ) :
    wallValue normal offset (segment a b t) =
      (1 - t) * wallValue normal offset a + t * wallValue normal offset b := by
  unfold wallValue segment
  simp only [mul_add, Finset.sum_add_distrib]
  have ha : (∑ i, normal i * ((1 - t) * a i)) = (1 - t) * ∑ i, normal i * a i := by
    rw [Finset.mul_sum]
    apply Finset.sum_congr rfl
    intro i _
    ring
  have hb : (∑ i, normal i * (t * b i)) = t * ∑ i, normal i * b i := by
    rw [Finset.mul_sum]
    apply Finset.sum_congr rfl
    intro i _
    ring
  rw [ha, hb]
  ring

theorem wall_hit_range (normal : Fin 3 → ℝ) (offset : ℝ) (a b : Fin 3 → ℝ)
    (ha : 0 < wallValue normal offset a) (hb : wallValue normal offset b < 0) :
    0 < wallHit normal offset a b ∧ wallHit normal offset a b < 1 := by
  have hd : 0 < wallValue normal offset a - wallValue normal offset b := by linarith
  unfold wallHit
  constructor
  · exact div_pos ha hd
  · apply (div_lt_one hd).2
    linarith

theorem wall_hit_on_surface (normal : Fin 3 → ℝ) (offset : ℝ) (a b : Fin 3 → ℝ)
    (ha : 0 < wallValue normal offset a) (hb : wallValue normal offset b < 0) :
    wallValue normal offset (segment a b (wallHit normal offset a b)) = 0 := by
  have hd : wallValue normal offset a - wallValue normal offset b ≠ 0 := by linarith
  rw [wall_segment_affine]
  unfold wallHit
  field_simp [hd]
  ring

theorem wall_first_hit (normal : Fin 3 → ℝ) (offset : ℝ) (a b : Fin 3 → ℝ)
    (ha : 0 < wallValue normal offset a) (hb : wallValue normal offset b < 0)
    (t : ℝ) (_ht0 : 0 ≤ t) (ht : t < wallHit normal offset a b) :
    0 < wallValue normal offset (segment a b t) := by
  have hd : 0 < wallValue normal offset a - wallValue normal offset b := by linarith
  have hmul : t * (wallValue normal offset a - wallValue normal offset b) <
      wallValue normal offset a := (lt_div_iff₀ hd).1 ht
  rw [wall_segment_affine]
  nlinarith

theorem clipped_segment_in_halfspace (normal : Fin 3 → ℝ) (offset : ℝ)
    (a b : Fin 3 → ℝ) (ha : 0 < wallValue normal offset a)
    (hb : wallValue normal offset b < 0) (s : ℝ) (_hs0 : 0 ≤ s) (hs1 : s ≤ 1) :
    0 ≤ wallValue normal offset (segment a (segment a b (wallHit normal offset a b)) s) := by
  rw [wall_segment_affine, wall_hit_on_surface normal offset a b ha hb]
  have h : 0 ≤ (1 - s) * wallValue normal offset a :=
    mul_nonneg (by linarith) (le_of_lt ha)
  simpa using h

def strain {n m : ℕ} (E : Fin m → Fin n → ℝ) (v : Fin n → ℝ) (e : Fin m) : ℝ :=
  ∑ i, E e i * v i

def viscousOperator {n m : ℕ} (E : Fin m → Fin n → ℝ) (mu : Fin m → ℝ)
    (v : Fin n → ℝ) (i : Fin n) : ℝ :=
  ∑ e, E e i * mu e * strain E v e

def dissipation {n m : ℕ} (E : Fin m → Fin n → ℝ) (mu : Fin m → ℝ)
    (v : Fin n → ℝ) : ℝ :=
  ∑ e, mu e * (strain E v e)^2

theorem viscous_work {n m : ℕ} (E : Fin m → Fin n → ℝ)
    (mu : Fin m → ℝ) (v : Fin n → ℝ) :
    (∑ i, viscousOperator E mu v i * v i) = dissipation E mu v := by
  unfold viscousOperator dissipation
  simp only [Finset.sum_mul]
  rw [Finset.sum_comm]
  apply Finset.sum_congr rfl
  intro e _
  calc
    (∑ i, E e i * mu e * strain E v e * v i) =
        ∑ i, (mu e * strain E v e) * (E e i * v i) := by
      apply Finset.sum_congr rfl
      intro i _
      ring
    _ = (mu e * strain E v e) * strain E v e := by
      rw [← Finset.mul_sum]
      rfl
    _ = mu e * (strain E v e)^2 := by ring

theorem dissipation_nonnegative {n m : ℕ} (E : Fin m → Fin n → ℝ)
    (mu : Fin m → ℝ) (v : Fin n → ℝ) (hm : ∀ e, 0 ≤ mu e) :
    0 ≤ dissipation E mu v := by
  exact Finset.sum_nonneg (fun e _ => mul_nonneg (hm e) (sq_nonneg _))

theorem backward_euler_work {n m : ℕ} (E : Fin m → Fin n → ℝ)
    (mu : Fin m → ℝ) (mass u v : Fin n → ℝ) (dt : ℝ)
    (step : ∀ i, mass i * v i + dt * viscousOperator E mu v i = mass i * u i) :
    (∑ i, mass i * (u i - v i) * v i) = dt * dissipation E mu v := by
  calc
    (∑ i, mass i * (u i - v i) * v i) =
        dt * ∑ i, viscousOperator E mu v i * v i := by
      rw [Finset.mul_sum]
      apply Finset.sum_congr rfl
      intro i _
      have hi : mass i * (u i - v i) = dt * viscousOperator E mu v i := by
        nlinarith [step i]
      rw [hi]
      ring
    _ = dt * dissipation E mu v := by rw [viscous_work]

theorem backward_euler_energy_nonincrease {n m : ℕ} (E : Fin m → Fin n → ℝ)
    (mu : Fin m → ℝ) (mass u v : Fin n → ℝ) (dt : ℝ)
    (hm : ∀ i, 0 ≤ mass i) (ht : 0 ≤ dt) (hmu : ∀ e, 0 ≤ mu e)
    (step : ∀ i, mass i * v i + dt * viscousOperator E mu v i = mass i * u i) :
    (∑ i, mass i * v i ^ 2) ≤ ∑ i, mass i * u i ^ 2 := by
  exact Rheon.Physics.implicit_energy_nonincrease mass u v dt (dissipation E mu v)
    hm ht (dissipation_nonnegative E mu v hmu) (backward_euler_work E mu mass u v dt step)

end
end Rheon.BoundedPhysics

Download Lean source

import Rheon.NoSlip

/-! Exact-real algebra for uniform tangential forcing in the fixed column.
Mass is rho times liquid dual volume; constraints are compatible and fixed.
Step, internal-force cancellation and assembled-power hypotheses are explicit.
These statements do not prove Rust/IEEE bounds, stencil assembly, the spectral
reference, transient convergence, mesh fluid coupling or capillarity. -/
namespace Rheon.ColumnForcing
noncomputable section
open scoped BigOperators

theorem units_agree (rho volume acceleration densityForce : ℝ)
    (units : densityForce = rho * acceleration) :
    (rho * volume) * acceleration = volume * densityForce := by rw [units]; ring

theorem reaction_includes_body_force (u wall mass dt viscous navier body : ℝ)
    (trace : u = wall) :
    u + (dt * (viscous + navier + body) + (-dt * (viscous + navier + body))) / mass = wall := by
  exact Rheon.NoSlip.compatible_reaction_preserves_trace u wall mass dt
    (viscous + navier + body) trace

theorem forced_momentum {n : ℕ} (mass delta internal body wall : Fin n → ℝ)
    (dt : ℝ)
    (step : ∀ i, mass i * delta i = dt * (internal i + body i) + wall i)
    (cancellation : ∑ i, internal i = 0) :
    (∑ i, mass i * delta i) = dt * (∑ i, body i) + ∑ i, wall i := by
  have h := Rheon.NoSlip.momentum_balance mass delta (fun i => internal i + body i) wall dt step
  simpa only [Finset.sum_add_distrib, cancellation, zero_add] using h

theorem forced_work {n : ℕ}
    (mass u delta force reaction wall : Fin n → ℝ)
    (dt bulk navierLoss navierPower bodyPower : ℝ)
    (step : ∀ i, mass i * delta i = dt * force i + reaction i)
    (power : (∑ i, force i * u i) = navierPower + bodyPower - bulk - navierLoss)
    (compatibleWork : ∀ i, reaction i * u i = reaction i * wall i) :
    (∑ i, mass i / 2 * ((u i + delta i)^2 - (u i)^2 - (delta i)^2)) +
      dt * (bulk + navierLoss - navierPower - bodyPower) -
      ∑ i, reaction i * wall i = 0 := by
  have h := Rheon.NoSlip.constrained_work_identity mass u delta force reaction wall
    dt bulk navierLoss (navierPower + bodyPower) step power compatibleWork
  convert h using 1; ring

def parabola (G mu H y : ℝ) : ℝ := G / (2 * mu) * y * (H - y)

theorem parabola_endpoints (G mu H : ℝ) :
    parabola G mu H 0 = 0 ∧ parabola G mu H H = 0 := by
  unfold parabola
  constructor <;> ring

theorem quadratic_stencil_balance (G mu H y h area : ℝ) (viscosity : mu ≠ 0)
    (spacing : h ≠ 0) :
    (mu * area / h) * (parabola G mu H (y-h) - 2 * parabola G mu H y +
      parabola G mu H (y+h)) + area * h * G = 0 := by
  unfold parabola
  field_simp
  ring

theorem recurrence_about_equilibrium (mass dt stiffness old next left right
    eq eqLeft eqRight body : ℝ)
    (step : mass * (next-old) = dt * (stiffness * (left-2*old+right) + body))
    (steady : stiffness * (eqLeft-2*eq+eqRight) + body = 0) :
    mass * (next-eq) = mass * (old-eq) +
      dt * stiffness * ((left-eqLeft)-2*(old-eq)+(right-eqRight)) := by
  nlinarith [mul_eq_zero_of_right dt steady]

end
end Rheon.ColumnForcing

Download Lean source

import Rheon.Discrete
import Rheon.Transport
import Mathlib.Algebra.Order.BigOperators.Ring.Finset
import Mathlib.Tactic.FieldSimp

/-! Exact-real finite-volume Euler contracts. Shared face fluxes are single values
with balanced incidence. Positivity is conditional on admitted states and explicit
convex/CFL coefficients. No Riemann solver, reconstruction, mesh assembly, IEEE
arithmetic, entropy inequality or continuum convergence is verified here. -/
namespace Rheon.Compressible
noncomputable section
open scoped BigOperators

def faceUpdate {n f k : ℕ} (B : Fin n → Fin f → ℝ)
    (amount : Fin n → Fin k → ℝ) (flux : Fin f → Fin k → ℝ)
    (dt : ℝ) (i : Fin n) (c : Fin k) : ℝ :=
  amount i c + dt * Rheon.Discrete.incidence B (fun e => flux e c) i

/-- One shared flux per face conserves each extensive Euler component. -/
theorem shared_face_conservation {n f k : ℕ} (B : Fin n → Fin f → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0)
    (amount : Fin n → Fin k → ℝ) (flux : Fin f → Fin k → ℝ)
    (dt : ℝ) (c : Fin k) :
    (∑ i, faceUpdate B amount flux dt i c) = ∑ i, amount i c := by
  exact Rheon.Discrete.conservative_update B balanced (fun i => amount i c)
    (fun e => flux e c) dt

def kineticDensity {d : ℕ} (rho : ℝ) (momentum : Fin d → ℝ) : ℝ :=
  ∑ j, (momentum j)^2 / (2 * rho)

def internalDensity {d : ℕ} (rho : ℝ) (momentum : Fin d → ℝ) (energy : ℝ) : ℝ :=
  energy - kineticDensity rho momentum

def pressure {d : ℕ} (gamma rho : ℝ) (momentum : Fin d → ℝ) (energy : ℝ) : ℝ :=
  (gamma - 1) * internalDensity rho momentum energy

def Admissible {d : ℕ} (rho : ℝ) (momentum : Fin d → ℝ) (energy : ℝ) : Prop :=
  0 < rho ∧ 0 < internalDensity rho momentum energy

theorem ideal_gas_pressure_positive {d : ℕ} (gamma rho : ℝ)
    (momentum : Fin d → ℝ) (energy : ℝ) (hg : 1 < gamma)
    (state : Admissible rho momentum energy) :
    0 < pressure gamma rho momentum energy := by
  exact mul_pos (sub_pos.mpr hg) state.2

theorem admissible_iff_pressure_positive {d : ℕ} (gamma rho : ℝ)
    (momentum : Fin d → ℝ) (energy : ℝ) (hg : 1 < gamma) :
    Admissible rho momentum energy ↔ 0 < rho ∧ 0 < pressure gamma rho momentum energy := by
  constructor
  · intro h
    exact ⟨h.1, ideal_gas_pressure_positive gamma rho momentum energy hg h⟩
  · rintro ⟨hr, hp⟩
    exact ⟨hr, (mul_pos_iff_of_pos_left (sub_pos.mpr hg)).mp hp⟩

/-- Positive endpoint values stay strictly positive, including t=0 and t=1. -/
theorem blend_positive (t a b : ℝ) (ht0 : 0 ≤ t) (ht1 : t ≤ 1)
    (ha : 0 < a) (hb : 0 < b) : 0 < Rheon.Transport.blend t a b := by
  unfold Rheon.Transport.blend
  by_cases h : t = 1
  · simp [h, hb]
  · have ht : t < 1 := lt_of_le_of_ne ht1 h
    have hp : 0 < (1 - t) * a := mul_pos (sub_pos.mpr ht) ha
    have hn : 0 ≤ t * b := mul_nonneg ht0 (le_of_lt hb)
    linarith

/-- The exact density-weighted Cauchy remainder; valid for signed momentum. -/
theorem kinetic_blend_remainder (t r s m n : ℝ) (hr : 0 < r) (hs : 0 < s)
    (ht0 : 0 ≤ t) (ht1 : t ≤ 1) :
    (1-t) * (m^2 / (2*r)) + t * (n^2 / (2*s)) -
      (Rheon.Transport.blend t m n)^2 / (2 * Rheon.Transport.blend t r s) =
    t * (1-t) * (s*m-r*n)^2 / (2*r*s*Rheon.Transport.blend t r s) := by
  have hb := blend_positive t r s ht0 ht1 hr hs
  unfold Rheon.Transport.blend at hb ⊢
  field_simp [ne_of_gt hr, ne_of_gt hs, ne_of_gt hb]
  ring

/-- Convexity is proved from the nonnegative square remainder, not assumed. -/
theorem kinetic_density_convex {d : ℕ} (t r s : ℝ) (m n : Fin d → ℝ)
    (hr : 0 < r) (hs : 0 < s) (ht0 : 0 ≤ t) (ht1 : t ≤ 1) :
    kineticDensity (Rheon.Transport.blend t r s) (fun j => Rheon.Transport.blend t (m j) (n j)) ≤
      Rheon.Transport.blend t (kineticDensity r m) (kineticDensity s n) := by
  have hb := blend_positive t r s ht0 ht1 hr hs
  unfold kineticDensity Rheon.Transport.blend
  rw [Finset.mul_sum, Finset.mul_sum, ← Finset.sum_add_distrib]
  apply Finset.sum_le_sum
  intro j _
  have hi := kinetic_blend_remainder t r s (m j) (n j) hr hs ht0 ht1
  have hn : 0 ≤ t * (1-t) * (s*m j-r*n j)^2 /
      (2*r*s*Rheon.Transport.blend t r s) :=
    div_nonneg (mul_nonneg (mul_nonneg ht0 (sub_nonneg.mpr ht1)) (sq_nonneg _))
      (le_of_lt (mul_pos (mul_pos (mul_pos (by norm_num) hr) hs) hb))
  unfold Rheon.Transport.blend at hi hn
  linarith

/-- The ideal-gas admissible set is convex in conservative variables. -/
theorem admissible_blend {d : ℕ} (t r s E F : ℝ) (m n : Fin d → ℝ)
    (ht0 : 0 ≤ t) (ht1 : t ≤ 1)
    (left : Admissible r m E) (right : Admissible s n F) :
    Admissible (Rheon.Transport.blend t r s)
      (fun j => Rheon.Transport.blend t (m j) (n j)) (Rheon.Transport.blend t E F) := by
  refine ⟨blend_positive t r s ht0 ht1 left.1 right.1, ?_⟩
  have hk := kinetic_density_convex t r s m n left.1 right.1 ht0 ht1
  have hp := blend_positive t (internalDensity r m E) (internalDensity s n F)
    ht0 ht1 left.2 right.2
  unfold internalDensity Rheon.Transport.blend at *
  linarith

def normalVelocity {d : ℕ} (rho : ℝ) (momentum normal : Fin d → ℝ) : ℝ :=
  (∑ j, momentum j * normal j) / rho

def eulerMassFlux (rho velocity : ℝ) : ℝ := rho * velocity

def eulerMomentumFlux {d : ℕ} (momentum normal : Fin d → ℝ)
    (p velocity : ℝ) (j : Fin d) : ℝ := momentum j * velocity + p * normal j

def eulerEnergyFlux (energy p velocity : ℝ) : ℝ := (energy + p) * velocity

/-- Stationary slip means zero normal momentum. Pressure still transfers momentum. -/
theorem stationary_slip_wall_flux {d : ℕ} (rho energy p : ℝ)
    (momentum normal : Fin d → ℝ) (slip : (∑ j, momentum j * normal j) = 0) :
    eulerMassFlux rho (normalVelocity rho momentum normal) = 0 ∧
    eulerEnergyFlux energy p (normalVelocity rho momentum normal) = 0 ∧
    ∀ j, eulerMomentumFlux momentum normal p (normalVelocity rho momentum normal) j = p * normal j := by
  simp [normalVelocity, slip, eulerMassFlux, eulerEnergyFlux, eulerMomentumFlux]

/-- A donor-density update with explicit outflow coefficients. This definition
is a convex density contract, not a general Euler Riemann-flux positivity claim. -/
def cflDensity {f : ℕ} (rho : ℝ) (coefficient donor : Fin f → ℝ) : ℝ :=
  (1 - ∑ e, coefficient e) * rho + ∑ e, coefficient e * donor e

theorem cfl_density_nonnegative {f : ℕ} (rho : ℝ) (coefficient donor : Fin f → ℝ)
    (hr : 0 ≤ rho) (hc : ∀ e, 0 ≤ coefficient e)
    (hd : ∀ e, 0 ≤ donor e) (cfl : (∑ e, coefficient e) ≤ 1) :
    0 ≤ cflDensity rho coefficient donor := by
  exact add_nonneg (mul_nonneg (sub_nonneg.mpr cfl) hr)
    (Finset.sum_nonneg (fun e _ => mul_nonneg (hc e) (hd e)))

/-- Strict CFL leaves a positive contribution from the current density. -/
theorem cfl_density_positive {f : ℕ} (rho : ℝ) (coefficient donor : Fin f → ℝ)
    (hr : 0 < rho) (hc : ∀ e, 0 ≤ coefficient e)
    (hd : ∀ e, 0 ≤ donor e) (cfl : (∑ e, coefficient e) < 1) :
    0 < cflDensity rho coefficient donor := by
  exact add_pos_of_pos_of_nonneg (mul_pos (sub_pos.mpr cfl) hr)
    (Finset.sum_nonneg (fun e _ => mul_nonneg (hc e) (hd e)))

end
end Rheon.Compressible

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Tactic.Ring
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.NormNum
namespace Rheon.CouplingStudy
noncomputable section
-- Finite affine sample geometry, not arbitrary mesh reconstruction.
theorem affine_transverse_sample (b t : ℝ) (ht : t ≠ 0) : b*t/t=b := by field_simp

theorem affine_longitudinal_sample (a b s : ℝ) : a+b*s-s*b=a := by ring

theorem nonorthogonal_correction (a b s : ℝ) (h : 1+s*s ≠ 0) :
 (a+b*s)/(1+s*s)+a*(1-1/(1+s*s))-b*s/(1+s*s)=a := by field_simp; ring

theorem skew_face_correction (b s c : ℝ) : c+b*s/2-b*s/2=c := by ring

-- Explicit recurrence of the stated two-path reduced momentum model.
theorem conservative_corrector (u v : ℝ) :
 (u-(u+v)/2)+(v-(u+v)/2)=0 := by ring

theorem velocity_error_recurrence (u a : ℝ) :
 ((1-a)*u+a*(1-2*u)/4)-1/6=(1-3*a/2)*(u-1/6) := by ring

theorem pressure_error_recurrence (p a : ℝ) :
 (p+a*(1/2-p))-1/2=(1-a)*(p-1/2) := by ring

theorem bounded_velocity_factor (a : ℝ) (ha : 0 ≤ a) (hb : a ≤ 4/3) :
 -1 ≤ 1-3*a/2 ∧ 1-3*a/2 ≤ 1 := by constructor <;> linarith

theorem exact_network_solution :
 (2*(1/6:ℝ)-(-1/6)+1/2=1) ∧ (2*(-1/6:ℝ)-1/6+1/2=0) := by norm_num

theorem overrelaxed_velocity_factor : (1-3*(3/2:ℝ)/2)^2=25/16 := by norm_num
end
end Rheon.CouplingStudy

Download Lean source

import Lean.Elab.Tactic.Omega

/-!
Exact natural-number contracts for Rheon's CPU decomposition and feedback.
These statements do not formalize Rust, Rayon, atomics, machine overflow,
wall-clock progress, numerical error, or a performance optimum.
-/
namespace Rheon.CpuScheduling

/-- The production clamp applied to each proposed linear batch size. -/
def batchSize (proposed : Nat) : Nat := min 16384 (max 64 proposed)

theorem batch_size_bounds (proposed : Nat) :
    64 ≤ batchSize proposed ∧ batchSize proposed ≤ 16384 := by
  unfold batchSize
  omega

/-- Membership in a clipped half-open scratch chunk; all arithmetic is exact. -/
def InChunk (cells grain chunk index : Nat) : Prop :=
  index < cells ∧ chunk * grain ≤ index ∧ index < (chunk + 1) * grain

theorem chunk_covers (cells grain index : Nat) (hg : 0 < grain) (hi : index < cells) :
    InChunk cells grain (index / grain) index := by
  have hdiv := Nat.div_add_mod index grain
  rw [Nat.mul_comm grain (index / grain)] at hdiv
  have hmod := Nat.mod_lt index hg
  unfold InChunk
  rw [Nat.add_mul]
  simp only [Nat.one_mul]
  omega

theorem chunk_unique (cells grain index a b : Nat)
    (ha : InChunk cells grain a index) (hb : InChunk cells grain b index) : a = b := by
  rcases ha with ⟨_, halower, haupper⟩
  rcases hb with ⟨_, hblower, hbupper⟩
  by_cases hab : a < b
  · have h : (a + 1) * grain ≤ b * grain := Nat.mul_le_mul_right grain hab
    omega
  · by_cases hba : b < a
    · have h : (b + 1) * grain ≤ a * grain := Nat.mul_le_mul_right grain hba
      omega
    · omega

/-- The wave rank of the three tile coordinates. -/
def waveRank (x y z : Nat) : Nat := x + y + z

/-- A forward triangular read crossing one tile boundary. -/
def Precedes (ax ay az bx by_ bz : Nat) : Prop :=
  (ax + 1 = bx ∧ ay = by_ ∧ az = bz) ∨
  (ax = bx ∧ ay + 1 = by_ ∧ az = bz) ∨
  (ax = bx ∧ ay = by_ ∧ az + 1 = bz)

theorem predecessor_rank (ax ay az bx by_ bz : Nat)
    (h : Precedes ax ay az bx by_ bz) :
    waveRank ax ay az + 1 = waveRank bx by_ bz := by
  unfold Precedes waveRank at *
  rcases h with h | h | h <;> omega

theorem same_wave_has_no_dependency (ax ay az bx by_ bz : Nat)
    (h : waveRank ax ay az = waveRank bx by_ bz) :
    ¬ Precedes ax ay az bx by_ bz ∧ ¬ Precedes bx by_ bz ax ay az := by
  constructor <;> intro hp
  · have := predecessor_rank ax ay az bx by_ bz hp
    omega
  · have := predecessor_rank bx by_ bz ax ay az hp
    omega

/-- The last rank of a nonempty rectangular tile grid. -/
def lastWave (nx ny nz : Nat) : Nat := nx + ny + nz - 3

theorem tile_wave_in_bounds (nx ny nz x y z : Nat)
    (hx : x < nx) (hy : y < ny) (hz : z < nz) :
    waveRank x y z ≤ lastWave nx ny nz := by
  unfold waveRank lastWave
  omega

theorem reverse_wave_dependency (ax ay az bx by_ bz last : Nat)
    (h : Precedes ax ay az bx by_ bz) (hb : waveRank bx by_ bz ≤ last) :
    last - waveRank bx by_ bz < last - waveRank ax ay az := by
  have := predecessor_rank ax ay az bx by_ bz h
  omega

/-- Exact phase policy; `completed` counts successful, nonoverlapping calls
of one stage with a stable workload, worker count and preferred grain. -/
def trialCandidate (completed : Nat) : Option Nat :=
  if completed % 128 < 8 then some ((completed % 128) % 4) else none

theorem trial_candidate_bounded (completed candidate : Nat)
    (h : trialCandidate completed = some candidate) : candidate < 4 := by
  unfold trialCandidate at h
  split at h
  · simp only [Option.some.injEq] at h
    omega
  · contradiction

theorem trial_first_round (candidate : Nat) (h : candidate < 4) :
    trialCandidate candidate = some candidate := by
  unfold trialCandidate
  have hp : candidate % 128 = candidate := by omega
  have hc : candidate % 4 = candidate := by omega
  simp [hp, hc, show candidate < 8 by omega]

theorem trial_second_round (candidate : Nat) (h : candidate < 4) :
    trialCandidate (candidate + 4) = some candidate := by
  unfold trialCandidate
  have hp : (candidate + 4) % 128 = candidate + 4 := by omega
  have hc : (candidate + 4) % 4 = candidate := by omega
  simp [hp, hc, show candidate + 4 < 8 by omega]

theorem trial_periodic (completed : Nat) :
    trialCandidate (completed + 128) = trialCandidate completed := by
  have h : (completed + 128) % 128 = completed % 128 := by omega
  simp only [trialCandidate, h]

theorem no_trial_after_exploration (completed : Nat) (h : 8 ≤ completed % 128) :
    trialCandidate completed = none := by
  simp [trialCandidate, show ¬ completed % 128 < 8 by omega]

theorem candidate_revisited_within_period (completed candidate : Nat) (h : candidate < 4) :
    ∃ delay : Nat, delay < 128 ∧ trialCandidate (completed + delay) = some candidate := by
  refine ⟨(candidate + 128 - completed % 128) % 128, ?_, ?_⟩
  · omega
  · have hp : (completed + (candidate + 128 - completed % 128) % 128) % 128 = candidate := by
      omega
    unfold trialCandidate
    simp [hp, show candidate < 8 by omega, show candidate % 4 = candidate by omega]

/-- A successful completion increments the phase counter once. Failed,
overlapping, reset or shape-changing operations do not satisfy this premise. -/
def SuccessfulSerialTrace (counter : Nat → Nat) : Prop :=
  ∀ call, counter (call + 1) = counter call + 1

theorem successful_trace_count (counter : Nat → Nat) (h : SuccessfulSerialTrace counter)
    (call : Nat) : counter call = counter 0 + call := by
  induction call with
  | zero => simp
  | succ call ih =>
      have hs := h call
      omega

theorem successful_trace_explores (counter : Nat → Nat) (h : SuccessfulSerialTrace counter)
    (start candidate : Nat) (hc : candidate < 4) :
    ∃ delay : Nat, delay < 128 ∧ trialCandidate (counter (start + delay)) = some candidate := by
  obtain ⟨delay, hd, ht⟩ := candidate_revisited_within_period (counter start) candidate hc
  refine ⟨delay, hd, ?_⟩
  have hs := successful_trace_count counter h start
  have he := successful_trace_count counter h (start + delay)
  have hcounters : counter (start + delay) = counter start + delay := by omega
  rw [hcounters]
  exact ht

end Rheon.CpuScheduling

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Algebra.BigOperators.Ring.Finset
import Mathlib.Algebra.Order.BigOperators.Group.Finset
import Mathlib.Data.Fintype.Fin
import Mathlib.Tactic.Ring

/-! Exact finite-dimensional contracts. B is a cells-by-faces incidence matrix
with +1 at a face head and -1 at its tail. Physical divergence is -B u.
No theorem in this module is a floating-point or continuum-flow theorem. -/
namespace Rheon.Discrete
noncomputable section

def gradient {n m : ℕ} (B : Fin n → Fin m → ℝ) (p : Fin n → ℝ) (e : Fin m) : ℝ :=
  ∑ i, B i e * p i

def incidence {n m : ℕ} (B : Fin n → Fin m → ℝ) (u : Fin m → ℝ) (i : Fin n) : ℝ :=
  ∑ e, B i e * u e

def laplace {n m : ℕ} (B : Fin n → Fin m → ℝ) (w : Fin m → ℝ)
    (p : Fin n → ℝ) (i : Fin n) : ℝ :=
  incidence B (fun e => w e * gradient B p e) i

def corrected {n m : ℕ} (B : Fin n → Fin m → ℝ) (w u : Fin m → ℝ)
    (p : Fin n → ℝ) (dt : ℝ) (e : Fin m) : ℝ :=
  u e - dt * w e * gradient B p e

theorem adjoint_identity {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (u : Fin m → ℝ) (p : Fin n → ℝ) :
    (∑ i, p i * incidence B u i) = ∑ e, gradient B p e * u e := by
  simp only [incidence, gradient, Finset.mul_sum, Finset.sum_mul]
  rw [Finset.sum_comm]
  apply Finset.sum_congr rfl
  intro e he
  apply Finset.sum_congr rfl
  intro i hi
  ring

theorem pressure_energy {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (w : Fin m → ℝ) (p : Fin n → ℝ) :
    (∑ i, p i * laplace B w p i) = ∑ e, w e * (gradient B p e)^2 := by
  unfold laplace
  rw [adjoint_identity]
  apply Finset.sum_congr rfl
  intro e he
  ring

theorem pressure_nonnegative {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (w : Fin m → ℝ) (p : Fin n → ℝ) (hw : ∀ e, 0 ≤ w e) :
    0 ≤ ∑ i, p i * laplace B w p i := by
  rw [pressure_energy]
  exact Finset.sum_nonneg (fun e _ => mul_nonneg (hw e) (sq_nonneg _))

theorem pressure_symmetric {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (w : Fin m → ℝ) (p q : Fin n → ℝ) :
    (∑ i, p i * laplace B w q i) = ∑ i, q i * laplace B w p i := by
  unfold laplace
  rw [adjoint_identity, adjoint_identity]
  apply Finset.sum_congr rfl
  intro e he
  ring

theorem constant_gradient_zero {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0) (c : ℝ) :
    gradient B (fun _ => c) = 0 := by
  funext e
  simp [gradient, ← Finset.sum_mul, balanced]

theorem constant_pressure_nullspace {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (w : Fin m → ℝ) (balanced : ∀ e, (∑ i, B i e) = 0) (c : ℝ) :
    laplace B w (fun _ => c) = 0 := by
  have h := constant_gradient_zero B balanced c
  funext i
  simp [laplace, incidence, h]

theorem internal_flux_conservation {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0) (u : Fin m → ℝ) :
    (∑ i, incidence B u i) = 0 := by
  unfold incidence
  rw [Finset.sum_comm]
  simp_rw [← Finset.sum_mul, balanced]
  simp

theorem conservative_update {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0) (q : Fin n → ℝ)
    (flux : Fin m → ℝ) (dt : ℝ) :
    (∑ i, (q i + dt * incidence B flux i)) = ∑ i, q i := by
  rw [Finset.sum_add_distrib, ← Finset.mul_sum, internal_flux_conservation B balanced]
  ring

theorem projection_residual {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (w u : Fin m → ℝ) (p b : Fin n → ℝ) (dt : ℝ)
    (rhs : ∀ i, incidence B u i = dt * b i) (i : Fin n) :
    incidence B (corrected B w u p dt) i = dt * (b i - laplace B w p i) := by
  unfold incidence corrected laplace
  simp only [mul_sub, Finset.sum_sub_distrib]
  have h : (∑ e, B i e * (dt * w e * gradient B p e)) =
      dt * ∑ e, B i e * (w e * gradient B p e) := by
    rw [Finset.mul_sum]
    apply Finset.sum_congr rfl
    intro e he
    ring
  rw [h]
  have hr := rhs i
  unfold incidence at hr
  rw [hr]
  unfold incidence
  ring

theorem exact_projection {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (w u : Fin m → ℝ) (p b : Fin n → ℝ) (dt : ℝ)
    (rhs : ∀ i, incidence B u i = dt * b i)
    (solved : ∀ i, laplace B w p i = b i) :
    incidence B (corrected B w u p dt) = 0 := by
  funext i
  rw [projection_residual B w u p b dt rhs i, solved]
  simp

end
end Rheon.Discrete

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Tactic.Ring
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Positivity

/-! Exact finite teaching contracts. Unit volumes, spacing and density are
specified by the lesson. No mesh reader, geometric quadrature, general
Rhie–Chow/SIMPLE convergence, PDE, Rust/JS or IEEE refinement is proved. -/
namespace Rheon.EngineeringFlux
noncomputable section

theorem shared_face_cancellation (flux component : ℝ) :
    flux * component + -(flux * component) = 0 := by ring

theorem two_cell_update_conserves (left right dt flux : ℝ) :
    (left - dt * flux) + (right + dt * flux) = left + right := by ring

theorem boundary_work (area velocity q pL pR : ℝ) :
    pL * (area * velocity + q) + pR * (-area * velocity - q) =
      -(pR - pL) * area * velocity + (pL - pR) * q := by ring

theorem balanced_flow (area velocity pL pR : ℝ) :
    -(pR - pL) * area * velocity + (pL - pR) * (-area * velocity) = 0 := by ring

theorem omitted_boundary_counterexample :
    (1 : ℝ) * (1/2 + 1/2) + 3 * (-1/2 - 1/2) = -2 ∧
    -(3 - 1) * (1/2 : ℝ) = -1 := by norm_num

theorem checkerboard_central_blind (a h : ℝ) :
    (-a - -a) / (2*h) = 0 ∧ (a-a) / (2*h) = 0 := by simp

theorem checkerboard_face_response (a r : ℝ) :
    -r * (-a-a) = 2*r*a ∧ -r * (a- -a) = -2*r*a := by constructor <;> ring

theorem checkerboard_divergence (a r : ℝ) :
    (2*r*a)-(-2*r*a) = 4*r*a := by ring

theorem periodic_pressure_energy (a b c d r : ℝ) :
    a*(r*(2*a-b-d)) + b*(r*(2*b-c-a)) +
      c*(r*(2*c-d-b)) + d*(r*(2*d-a-c)) =
      r*((b-a)^2+(c-b)^2+(d-c)^2+(a-d)^2) := by ring

def planarPower (mu a b c d : ℝ) : ℝ :=
    2*mu*a^2 + 2*mu*d^2 + mu*(b+c)^2

theorem planar_power_nonnegative (mu a b c d : ℝ) (hmu : 0 ≤ mu) :
    0 ≤ planarPower mu a b c d := by unfold planarPower; positivity

-- Face-centre quadrature for the specified centred unit cube reduces to
-- tau:gradient. Quadrature correctness for arbitrary geometry is not asserted.
theorem planar_boundary_power (mu a b c d : ℝ) :
    (2*mu*a)*a + (mu*(b+c))*b + (mu*(b+c))*c + (2*mu*d)*d =
      planarPower mu a b c d := by unfold planarPower; ring

theorem rigid_rotation_power (mu rate : ℝ) :
    planarPower mu 0 (-rate) rate 0 = 0 := by unfold planarPower; ring

theorem simple_shear_power (mu rate : ℝ) :
    planarPower mu 0 rate 0 0 = mu*rate^2 := by unfold planarPower; ring

theorem planar_extension_power (mu rate : ℝ) :
    planarPower mu rate 0 0 (-rate) = 4*mu*rate^2 := by unfold planarPower; ring

end
end Rheon.EngineeringFlux

Download Lean source

import Mathlib.Tactic.Ring

namespace Rheon.Indexing

def flatten (nx ny i j k : ℕ) : ℕ := i + nx * (j + ny * k)

theorem flatten_in_bounds (nx ny nz i j k : ℕ)
    (hi : i < nx) (hj : j < ny) (hk : k < nz) :
    flatten nx ny i j k < nx * ny * nz := by
  unfold flatten
  have hinner : j + ny * k < ny * nz := by
    calc
      j + ny * k < ny + ny * k := Nat.add_lt_add_right hj _
      _ = ny * (k + 1) := by ring
      _ ≤ ny * nz := Nat.mul_le_mul_left ny hk
  calc
    i + nx * (j + ny * k) < nx + nx * (j + ny * k) := Nat.add_lt_add_right hi _
    _ = nx * (j + ny * k + 1) := by ring
    _ ≤ nx * (ny * nz) := Nat.mul_le_mul_left nx hinner
    _ = nx * ny * nz := by ring

theorem staggered_count (n : ℕ) :
    (n+1)*n*n + n*(n+1)*n + n*n*(n+1) = 3*n*n*(n+1) := by ring

end Rheon.Indexing

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Data.Real.Sqrt
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.Ring

/-!
Exact-real contracts for the bounded material preview, not a verification of
`ParticleFluidSimulation`, its floating-point execution, neighbor enumeration,
time-step controller, or continuum/multiphase physics.

The pair coefficients below are the coordinate formulas in
`src/particle_fluid.rs`: reference densities `ri,rj`, estimated densities
`di,dj`, particle masses `mi,mj`, and the shared radial kernel factor. Their
symmetry is derived from those formulas, not assumed. The pressure statements
do not assert a compression-energy derivative across the 0.05 density floor.
Rust evaluates radial forces only for `1e-10*h ≤ distance < h`; the polynomial
below specifies that active-pair factor, not an extension of the force outside
its admitted support. Neighbor/support assembly is not proved here.
The viscosity energy bound is for one isolated pair without pressure, source,
gravity, cohesion, wall contact or changing coefficients. The actual controller
has additional conditions; these statements do not prove its global stability.

The nozzle profile/head and reservoir ledger describe the prescribed source.
They do not validate the approximate 0.43 disk-average factor, gas choking,
surface tension, temperature-dependent properties, or wall-pressure support.
-/
namespace Rheon.MaterialPreview
noncomputable section

def liquidPressure (c density reference : ℝ) : ℝ :=
  max (c ^ 2 * (density - reference)) 0

def gasPressure (density gasConstant temperature : ℝ) : ℝ :=
  density * gasConstant * temperature

theorem liquid_pressure_nonnegative (c density reference : ℝ) :
    0 ≤ liquidPressure c density reference := le_max_right _ _

theorem liquid_pressure_below_reference (c density reference : ℝ)
    (h : density ≤ reference) : liquidPressure c density reference = 0 := by
  apply max_eq_right
  exact mul_nonpos_of_nonneg_of_nonpos (sq_nonneg c) (sub_nonpos.mpr h)

theorem gas_reference_pressure (pressure gasConstant temperature : ℝ)
    (hr : gasConstant ≠ 0) (ht : temperature ≠ 0) :
    gasPressure (pressure / (gasConstant * temperature)) gasConstant temperature =
      pressure := by
  unfold gasPressure
  field_simp <;> ring

theorem reference_volume_contribution (reference volume neighborReference kernel : ℝ)
    (hn : neighborReference ≠ 0) :
    reference * ((neighborReference * volume) / neighborReference * kernel) =
      reference * volume * kernel := by
  field_simp <;> ring

def pressurePairCoefficient (pi pj ri rj di dj : ℝ) : ℝ :=
  pi * ri / (rj * di ^ 2) + pj * rj / (ri * dj ^ 2)

def pressurePairForce (mi mj pi pj ri rj di dj gradient displacement : ℝ) : ℝ :=
  -mi * mj * pressurePairCoefficient pi pj ri rj di dj * gradient * displacement

theorem pressure_pair_coefficient_symmetric (pi pj ri rj di dj : ℝ) :
    pressurePairCoefficient pi pj ri rj di dj =
      pressurePairCoefficient pj pi rj ri dj di := by
  unfold pressurePairCoefficient
  ring

theorem pressure_pair_action_reaction (mi mj pi pj ri rj di dj gradient delta : ℝ) :
    pressurePairForce mi mj pi pj ri rj di dj gradient delta +
      pressurePairForce mj mi pj pi rj ri dj di gradient (-delta) = 0 := by
  unfold pressurePairForce
  rw [pressure_pair_coefficient_symmetric pj pi rj ri dj di]
  ring

def poly6RadialGradient (coefficient h distance : ℝ) : ℝ :=
  -6 * coefficient * (h ^ 2 - distance ^ 2) ^ 2

theorem poly6_radial_gradient_nonpositive (coefficient h distance : ℝ)
    (hc : 0 ≤ coefficient) : poly6RadialGradient coefficient h distance ≤ 0 := by
  unfold poly6RadialGradient
  nlinarith [mul_nonneg hc (sq_nonneg (h ^ 2 - distance ^ 2))]

def viscosityPairCoefficient (mi mj mui muj di dj laplacian : ℝ) : ℝ :=
  mi * mj * ((mui + muj) / 2) / (di * dj) * laplacian

theorem viscosity_pair_coefficient_symmetric (mi mj mui muj di dj laplacian : ℝ) :
    viscosityPairCoefficient mi mj mui muj di dj laplacian =
      viscosityPairCoefficient mj mi muj mui dj di laplacian := by
  unfold viscosityPairCoefficient
  congr 1
  congr 1 <;> ring

theorem viscosity_pair_coefficient_nonnegative (mi mj mui muj di dj laplacian : ℝ)
    (hi : 0 ≤ mi) (hj : 0 ≤ mj) (hui : 0 ≤ mui) (huj : 0 ≤ muj)
    (hdi : 0 < di) (hdj : 0 < dj) (hl : 0 ≤ laplacian) :
    0 ≤ viscosityPairCoefficient mi mj mui muj di dj laplacian := by
  unfold viscosityPairCoefficient
  exact mul_nonneg (div_nonneg
    (mul_nonneg (mul_nonneg hi hj) (div_nonneg (add_nonneg hui huj) (by norm_num)))
    (le_of_lt (mul_pos hdi hdj))) hl

def viscosityForce (coefficient ui uj : ℝ) : ℝ := coefficient * (uj - ui)

theorem viscosity_pair_work (coefficient ui uj : ℝ) :
    ui * viscosityForce coefficient ui uj + uj * viscosityForce coefficient uj ui =
      -coefficient * (ui - uj) ^ 2 := by
  unfold viscosityForce
  ring

theorem viscosity_pair_work_nonpositive (coefficient ui uj : ℝ)
    (hc : 0 ≤ coefficient) :
    ui * viscosityForce coefficient ui uj + uj * viscosityForce coefficient uj ui ≤ 0 := by
  rw [viscosity_pair_work]
  have loss := mul_nonneg hc (sq_nonneg (ui - uj))
  nlinarith

def pairVelocity (mass dt coefficient ui uj : ℝ) : ℝ :=
  ui + dt * viscosityForce coefficient ui uj / mass

theorem viscosity_pair_momentum (mi mj dt coefficient ui uj : ℝ)
    (hi : mi ≠ 0) (hj : mj ≠ 0) :
    mi * pairVelocity mi dt coefficient ui uj +
      mj * pairVelocity mj dt coefficient uj ui = mi * ui + mj * uj := by
  unfold pairVelocity viscosityForce
  field_simp <;> ring

theorem viscosity_pair_energy_identity (mi mj dt coefficient ui uj : ℝ)
    (hi : mi ≠ 0) (hj : mj ≠ 0) :
    mi / 2 * (pairVelocity mi dt coefficient ui uj) ^ 2 +
      mj / 2 * (pairVelocity mj dt coefficient uj ui) ^ 2 -
      (mi / 2 * ui ^ 2 + mj / 2 * uj ^ 2) =
      -(dt * coefficient) *
        (1 - dt * coefficient / 2 * (1 / mi + 1 / mj)) * (ui - uj) ^ 2 := by
  unfold pairVelocity viscosityForce
  field_simp <;> ring

theorem viscosity_pair_energy_nonincrease (mi mj dt coefficient ui uj : ℝ)
    (hi : 0 < mi) (hj : 0 < mj) (hdt : 0 ≤ dt) (hc : 0 ≤ coefficient)
    (hstep : dt * coefficient * (1 / mi + 1 / mj) ≤ 2) :
    mi / 2 * (pairVelocity mi dt coefficient ui uj) ^ 2 +
      mj / 2 * (pairVelocity mj dt coefficient uj ui) ^ 2 ≤
      mi / 2 * ui ^ 2 + mj / 2 * uj ^ 2 := by
  have identity := viscosity_pair_energy_identity mi mj dt coefficient ui uj
    (ne_of_gt hi) (ne_of_gt hj)
  have factor : 0 ≤ 1 - dt * coefficient / 2 * (1 / mi + 1 / mj) := by
    nlinarith
  have loss := mul_nonneg (mul_nonneg (mul_nonneg hdt hc) factor) (sq_nonneg (ui - uj))
  nlinarith

def nozzleProfile (radialSquared verticalFraction : ℝ) : ℝ :=
  (1 - radialSquared) ^ 2 * (3 / 4 + verticalFraction / 4)

theorem nozzle_profile_bounds (radialSquared verticalFraction : ℝ)
    (hr0 : 0 ≤ radialSquared) (hr1 : radialSquared ≤ 1)
    (hv0 : -1 ≤ verticalFraction) (hv1 : verticalFraction ≤ 1) :
    0 ≤ nozzleProfile radialSquared verticalFraction ∧
      nozzleProfile radialSquared verticalFraction ≤ 1 := by
  have hs0 := sq_nonneg (1 - radialSquared)
  have hs1 : (1 - radialSquared) ^ 2 ≤ 1 := by nlinarith
  have ht0 : 0 ≤ (3 / 4 : ℝ) + verticalFraction / 4 := by linarith
  have ht1 : (3 / 4 : ℝ) + verticalFraction / 4 ≤ 1 := by linarith
  unfold nozzleProfile
  constructor
  · exact mul_nonneg hs0 ht0
  · nlinarith [mul_nonneg (sub_nonneg.mpr hs1) ht0]

def nozzleSpeed (head density profile : ℝ) : ℝ :=
  Real.sqrt (2 * head / density) * Real.sqrt profile

theorem nozzle_pressure_head_identity (head density profile : ℝ)
    (hh : 0 ≤ head) (hd : 0 < density) (hp : 0 ≤ profile) :
    density / 2 * (nozzleSpeed head density profile) ^ 2 = head * profile := by
  unfold nozzleSpeed
  rw [mul_pow, Real.sq_sqrt (div_nonneg (by linarith) (le_of_lt hd)), Real.sq_sqrt hp]
  field_simp <;> ring

theorem reservoir_mass_ledger (density credit flux dt parcelVolume count : ℝ) :
    count * (density * parcelVolume) +
      density * (credit + flux * dt - count * parcelVolume) =
      density * (credit + flux * dt) := by ring

end
end Rheon.MaterialPreview

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Ring
import Mathlib.Tactic.FieldSimp

/-! Exact amplification contracts for one scalar viscous mode.
`a` represents lambda * dt; energy is the unit cube's shear-mode energy for
rho=1. The PDE ansatz, trigonometric integral, exponential solution, sampled
3D rendering, Rust/JS execution and IEEE behavior are not formalized here. -/
namespace Rheon.ModalDecay
noncomputable section

def energy (amplitude : ℝ) : ℝ := amplitude ^ 2 / 4
def forward (a amplitude : ℝ) : ℝ := (1-a)*amplitude
def backward (a amplitude : ℝ) : ℝ := amplitude/(1+a)

theorem forward_energy_identity (a amplitude : ℝ) :
    energy (forward a amplitude) - energy amplitude =
      a*(a-2)*energy amplitude := by
  unfold energy forward
  ring

theorem forward_energy_nonincrease (a amplitude : ℝ)
    (ha : 0 ≤ a) (ha2 : a ≤ 2) :
    energy (forward a amplitude) ≤ energy amplitude := by
  have hprod : 0 ≤ a*(2-a) := mul_nonneg ha (by linarith)
  have h := mul_nonneg hprod (sq_nonneg amplitude)
  unfold energy forward
  nlinarith

theorem forward_sign_preserved (a amplitude : ℝ)
    (ha1 : a ≤ 1) (hamp : 0 ≤ amplitude) : 0 ≤ forward a amplitude := by
  exact mul_nonneg (by linarith) hamp

theorem forward_reverses (a amplitude : ℝ)
    (ha : 1 < a) (hamp : 0 < amplitude) : forward a amplitude < 0 := by
  exact mul_neg_of_neg_of_pos (by linarith) hamp

theorem backward_satisfies_update (a amplitude : ℝ) (ha : 0 ≤ a) :
    (1+a)*backward a amplitude = amplitude := by
  have hn : (1+a) ≠ 0 := by linarith
  unfold backward
  field_simp [hn]

theorem backward_energy_identity (a old next : ℝ)
    (update : (1+a)*next = old) :
    energy next - energy old = -(a*(2+a))*energy next := by
  rw [← update]
  unfold energy
  ring

theorem backward_energy_nonincrease (a old next : ℝ)
    (ha : 0 ≤ a) (update : (1+a)*next = old) : energy next ≤ energy old := by
  have hprod : 0 ≤ a*(2+a) := mul_nonneg ha (by linarith)
  have h := mul_nonneg hprod (sq_nonneg next)
  have identity := backward_energy_identity a old next update
  unfold energy at *
  nlinarith

theorem stability_does_not_preserve_sign :
    forward 2 1 = -1 ∧ energy (forward 2 1) = energy 1 := by
  norm_num [forward, energy]

end
end Rheon.ModalDecay

Download Lean source

import Rheon.WallFriction

/-! Exact-real identities for compatible, fixed endpoint constraints.
The reaction impulse cancels the unconstrained force impulse at constrained
nodes. These statements assume the step equation and assembled force power;
they do not prove stencil assembly, a Rust/IEEE bound, stability or convergence.
An incompatible initial trace or a changing prescribed speed is not modeled. -/
namespace Rheon.NoSlip
noncomputable section
open scoped BigOperators

theorem compatible_reaction_preserves_trace (u wall mass dt force : ℝ)
    (trace : u = wall) :
    u + (dt * force + (-dt * force)) / mass = wall := by
  rw [trace]
  ring

theorem relative_reaction_work_zero (u wall impulse : ℝ) (trace : u = wall) :
    impulse * (u - wall) = 0 := by rw [trace]; ring

theorem momentum_balance {n : ℕ} (mass delta force impulse : Fin n → ℝ)
    (dt : ℝ)
    (step : ∀ i, mass i * delta i = dt * force i + impulse i) :
    (∑ i, mass i * delta i) = dt * (∑ i, force i) + ∑ i, impulse i := by
  calc
    _ = ∑ i, (dt * force i + impulse i) := by
      apply Finset.sum_congr rfl
      intro i _
      exact step i
    _ = dt * (∑ i, force i) + ∑ i, impulse i := by
      rw [Finset.sum_add_distrib, Finset.mul_sum]

theorem constrained_work_identity {n : ℕ}
    (mass u delta force impulse wall : Fin n → ℝ)
    (dt bulk navierDissipation navierPower : ℝ)
    (step : ∀ i, mass i * delta i = dt * force i + impulse i)
    (power : (∑ i, force i * u i) = navierPower - bulk - navierDissipation)
    (compatibleWork : ∀ i, impulse i * u i = impulse i * wall i) :
    (∑ i, mass i / 2 * ((u i + delta i)^2 - (u i)^2 - (delta i)^2)) +
      dt * (bulk + navierDissipation - navierPower) - ∑ i, impulse i * wall i = 0 := by
  have hp : (∑ i, mass i * delta i * u i) =
      dt * (navierPower - bulk - navierDissipation) + ∑ i, impulse i * wall i := by
    calc
      _ = ∑ i, (dt * (force i * u i) + impulse i * wall i) := by
        apply Finset.sum_congr rfl
        intro i _
        rw [step]
        calc
          (dt * force i + impulse i) * u i =
              dt * (force i * u i) + impulse i * u i := by ring
          _ = dt * (force i * u i) + impulse i * wall i := by rw [compatibleWork]
      _ = dt * (∑ i, force i * u i) + ∑ i, impulse i * wall i := by
        rw [Finset.sum_add_distrib, Finset.mul_sum]
      _ = dt * (navierPower - bulk - navierDissipation) + ∑ i, impulse i * wall i := by rw [power]
  calc
    _ = (∑ i, mass i * delta i * u i) +
        dt * (bulk + navierDissipation - navierPower) - ∑ i, impulse i * wall i := by
      congr 2
      apply Finset.sum_congr rfl
      intro i _
      ring
    _ = 0 := by rw [hp]; ring

end
end Rheon.NoSlip

Download Lean source

import Rheon.Discrete
import Rheon.Physics

/-! Exact-real reference contracts for stationary obstacle pressure and reduced
extruded shear. `edgeIncidence` is +1 at the head and -1 at the tail; outward
volume flux is its negative applied to area times speed. Pressure inertia is
the declared face-area mass rho*A*d, not exact clipped dual volume. Shear uses
rho*V layer masses, shared-area/centroid-distance conductances and stationary
wall diagonals. The scalar shear interpretation is restricted to u=(u(y),0,0)
with constant Newtonian coefficients and invariant/periodic tangential axes.

Geometry admission and coefficient assembly are external obligations. These
statements do not certify Rust, IEEE arithmetic, PCG/Thomas convergence,
existence/uniqueness, complete nullspace characterization, general 3D viscous
traction, continuum convergence or moving walls. -/
namespace Rheon.ObstacleOperators
noncomputable section
open scoped BigOperators
open Rheon.Discrete

/-- One shared oriented face; tail and head are retained wet-cell indices. -/
def edgeIncidence {n m : ℕ} (tail head : Fin m → Fin n)
    (i : Fin n) (e : Fin m) : ℝ :=
  (if i = head e then 1 else 0) - (if i = tail e then 1 else 0)

theorem edge_gradient {n m : ℕ} (tail head : Fin m → Fin n)
    (p : Fin n → ℝ) (e : Fin m) :
    gradient (edgeIncidence tail head) p e = p (head e) - p (tail e) := by
  simp [gradient, edgeIncidence, sub_mul, Finset.sum_sub_distrib, ite_mul]

theorem edge_balanced {n m : ℕ} (tail head : Fin m → Fin n) (e : Fin m) :
    (∑ i, edgeIncidence tail head i e) = 0 := by
  simp [edgeIncidence, Finset.sum_sub_distrib]

def outwardFlux {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (area u : Fin m → ℝ) (i : Fin n) : ℝ :=
  -incidence B (fun e => area e * u e) i

def pressureWeight {m : ℕ} (rho : ℝ) (area distance : Fin m → ℝ)
    (e : Fin m) : ℝ := area e / (rho * distance e)

def faceMass {m : ℕ} (rho : ℝ) (area distance : Fin m → ℝ)
    (e : Fin m) : ℝ := rho * area e * distance e

def pressureCorrected {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (rho dt : ℝ) (distance u : Fin m → ℝ) (p : Fin n → ℝ) (e : Fin m) : ℝ :=
  u e - dt * gradient B p e / (rho * distance e)

def pressureRhs {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (area u : Fin m → ℝ) (dt : ℝ) (i : Fin n) : ℝ :=
  -outwardFlux B area u i / dt

def pressureResidual {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (rho dt : ℝ) (area distance u : Fin m → ℝ) (p : Fin n → ℝ) (i : Fin n) : ℝ :=
  pressureRhs B area u dt i - laplace B (pressureWeight rho area distance) p i

def kineticEnergy {n : ℕ} (mass u : Fin n → ℝ) : ℝ :=
  ∑ i, mass i / 2 * u i ^ 2

/-- This algebraic identity makes the implicit/new-time work convention explicit. -/
theorem kinetic_increment {n : ℕ} (mass u v : Fin n → ℝ) :
    kineticEnergy mass v - kineticEnergy mass u +
      kineticEnergy mass (fun i => v i - u i) =
      ∑ i, mass i * (v i - u i) * v i := by
  unfold kineticEnergy
  rw [← Finset.sum_sub_distrib, ← Finset.sum_add_distrib]
  apply Finset.sum_congr rfl
  intro i _
  ring

theorem face_mass_coefficient {m : ℕ} (rho : ℝ) (area distance : Fin m → ℝ)
    (hr : rho ≠ 0) (hd : ∀ e, distance e ≠ 0) (e : Fin m) :
    faceMass rho area distance e / (rho * distance e) = area e := by
  unfold faceMass
  field_simp [hr, hd e]
  ring

theorem pressure_corrected_flux {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (rho dt : ℝ) (area distance u : Fin m → ℝ) (p : Fin n → ℝ) :
    (fun e => area e * pressureCorrected B rho dt distance u p e) =
      corrected B (pressureWeight rho area distance) (fun e => area e * u e) p dt := by
  funext e
  unfold pressureCorrected corrected pressureWeight
  ring

/-- r=b-Lp, b=-Q_old/dt imply Qnew=-dt*r, including all gauge rows. -/
theorem pressure_flux_residual {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (rho dt : ℝ) (area distance u : Fin m → ℝ) (p : Fin n → ℝ)
    (ht : dt ≠ 0) (i : Fin n) :
    outwardFlux B area (pressureCorrected B rho dt distance u p) i =
      -dt * pressureResidual B rho dt area distance u p i := by
  have rhs : ∀ i, incidence B (fun e => area e * u e) i =
      dt * pressureRhs B area u dt i := by
    intro i
    unfold pressureRhs outwardFlux
    simp only [neg_neg]
    field_simp [ht]
  unfold outwardFlux
  rw [pressure_corrected_flux, projection_residual B _ _ _ _ dt rhs i]
  unfold pressureResidual
  ring

/-- A positive wet volume is essential when converting an integrated residual. -/
theorem pressure_divergence_residual {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (rho dt : ℝ) (area distance u : Fin m → ℝ) (p volume : Fin n → ℝ)
    (ht : 0 < dt) (hv : ∀ i, 0 < volume i) (i : Fin n) :
    outwardFlux B area (pressureCorrected B rho dt distance u p) i / volume i =
      -dt * pressureResidual B rho dt area distance u p i / volume i ∧
    |outwardFlux B area (pressureCorrected B rho dt distance u p) i / volume i| =
      dt * |pressureResidual B rho dt area distance u p i| / volume i := by
  rw [pressure_flux_residual B rho dt area distance u p (ne_of_gt ht) i]
  constructor
  · rfl
  · rw [abs_div, abs_mul, abs_neg, abs_of_pos ht, abs_of_pos (hv i)]

/-- Compatibility is local to each supplied connected label, not just global. -/
theorem component_rhs_compatible {n m : ℕ} (tail head : Fin m → Fin n)
    (label : Fin n → ℕ) (connected : ∀ e, label (tail e) = label (head e))
    (area u : Fin m → ℝ) (dt : ℝ) (c : ℕ) :
    (∑ i, (if label i = c then (1 : ℝ) else 0) *
      pressureRhs (edgeIncidence tail head) area u dt i) = 0 := by
  let potential : Fin n → ℝ := fun i => if label i = c then 1 else 0
  have hgrad : ∀ e, gradient (edgeIncidence tail head) potential e = 0 := by
    intro e
    rw [edge_gradient]
    simp only [potential, connected e, sub_self]
  have h := adjoint_identity (edgeIncidence tail head) (fun e => area e * u e) potential
  simp only [hgrad, zero_mul, Finset.sum_const_zero] at h
  unfold pressureRhs outwardFlux
  simp only [neg_neg, div_eq_mul_inv, ← mul_assoc, ← Finset.sum_mul]
  change (∑ i, potential i * incidence (edgeIncidence tail head) (fun e => area e * u e) i) * dt⁻¹ = 0
  rw [h, zero_mul]

/-- Derived from the coordinate correction and rho*A*d inertia, without a solve assumption. -/
theorem pressure_correction_energy {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (rho dt : ℝ) (area distance u : Fin m → ℝ) (p : Fin n → ℝ)
    (hr : rho ≠ 0) (hd : ∀ e, distance e ≠ 0) :
    kineticEnergy (faceMass rho area distance) (pressureCorrected B rho dt distance u p) -
      kineticEnergy (faceMass rho area distance) u +
      kineticEnergy (faceMass rho area distance)
        (fun e => pressureCorrected B rho dt distance u p e - u e) =
      dt * ∑ i, p i * outwardFlux B area (pressureCorrected B rho dt distance u p) i := by
  rw [kinetic_increment]
  have local_step : ∀ e, faceMass rho area distance e *
      (pressureCorrected B rho dt distance u p e - u e) =
      -dt * gradient B p e * area e := by
    intro e
    unfold pressureCorrected faceMass
    field_simp [hr, hd e]
    ring
  simp_rw [local_step]
  unfold outwardFlux
  simp only [mul_neg, Finset.sum_neg_distrib]
  rw [adjoint_identity]
  rw [Finset.mul_sum, ← Finset.sum_neg_distrib]
  apply Finset.sum_congr rfl
  intro e _
  ring

theorem exact_pressure_energy_nonincrease {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (rho dt : ℝ) (area distance u : Fin m → ℝ) (p : Fin n → ℝ)
    (hr : 0 < rho) (ha : ∀ e, 0 < area e) (hd : ∀ e, 0 < distance e)
    (ht : dt ≠ 0)
    (solved : ∀ i, pressureResidual B rho dt area distance u p i = 0) :
    kineticEnergy (faceMass rho area distance) (pressureCorrected B rho dt distance u p) ≤
      kineticEnergy (faceMass rho area distance) u := by
  have he := pressure_correction_energy B rho dt area distance u p
    (ne_of_gt hr) (fun e => ne_of_gt (hd e))
  simp_rw [pressure_flux_residual B rho dt area distance u p ht, solved, mul_zero] at he
  simp only [Finset.sum_const_zero, mul_zero] at he
  have inc : 0 ≤ kineticEnergy (faceMass rho area distance)
      (fun e => pressureCorrected B rho dt distance u p e - u e) := by
    apply Finset.sum_nonneg
    intro e _
    exact mul_nonneg (div_nonneg (le_of_lt (mul_pos (mul_pos hr (ha e)) (hd e)))
      (by norm_num)) (sq_nonneg _)
  linarith

/-- mu*A/d for shared fluid interfaces or stationary no-slip wall half-layers. -/
def shearConductance (mu area distance : ℝ) : ℝ := mu * area / distance

/-- Eliminated Navier wall trace; beta=0 is exactly free slip. -/
def navierConductance (mu area beta distance : ℝ) : ℝ :=
  area * mu * beta / (mu + beta * distance)

def navierTrace (mu beta distance speed : ℝ) : ℝ :=
  mu * speed / (mu + beta * distance)

theorem shear_conductance_positive (mu area distance : ℝ)
    (hm : 0 < mu) (ha : 0 < area) (hd : 0 < distance) :
    0 < shearConductance mu area distance := div_pos (mul_pos hm ha) hd

theorem navier_conductance_nonnegative (mu area beta distance : ℝ)
    (hm : 0 < mu) (ha : 0 ≤ area) (hb : 0 ≤ beta) (hd : 0 ≤ distance) :
    0 ≤ navierConductance mu area beta distance := by
  exact div_nonneg (mul_nonneg (mul_nonneg ha (le_of_lt hm)) hb)
    (le_of_lt (add_pos_of_pos_of_nonneg hm (mul_nonneg hb hd)))

theorem navier_free_slip (mu area distance : ℝ) :
    navierConductance mu area 0 distance = 0 := by
  simp [navierConductance]

/-- The effective force equals both fluid half-layer traction and wall friction. -/
theorem navier_series_traction (mu area beta distance speed : ℝ)
    (hm : 0 < mu) (hb : 0 ≤ beta) (hd : 0 < distance) :
    shearConductance mu area distance * (speed - navierTrace mu beta distance speed) =
      navierConductance mu area beta distance * speed ∧
    area * beta * navierTrace mu beta distance speed =
      navierConductance mu area beta distance * speed := by
  have hden : mu + beta * distance ≠ 0 :=
    ne_of_gt (add_pos_of_pos_of_nonneg hm (mul_nonneg hb (le_of_lt hd)))
  have hd0 : distance ≠ 0 := ne_of_gt hd
  unfold shearConductance navierTrace navierConductance
  constructor <;> field_simp [hden, hd0] <;> ring

/-- Boundary loss includes fluid half-layer strain and stationary-wall friction. -/
theorem navier_dissipation_split (mu area beta distance speed : ℝ)
    (hm : 0 < mu) (hb : 0 ≤ beta) (hd : 0 < distance) :
    shearConductance mu area distance * (speed - navierTrace mu beta distance speed)^2 +
      area * beta * (navierTrace mu beta distance speed)^2 =
      navierConductance mu area beta distance * speed^2 := by
  have hden : mu + beta * distance ≠ 0 :=
    ne_of_gt (add_pos_of_pos_of_nonneg hm (mul_nonneg hb (le_of_lt hd)))
  have hd0 : distance ≠ 0 := ne_of_gt hd
  unfold shearConductance navierTrace navierConductance
  field_simp [hden, hd0]
  ring

/-- Two wall contributions add even when both touch the same singleton layer. -/
def stationaryWallDiagonal {n : ℕ} (lower upper : Fin n) (kl ku : ℝ) (i : Fin n) : ℝ :=
  (if i = lower then kl else 0) + (if i = upper then ku else 0)

theorem stationary_wall_nonnegative {n : ℕ} (lower upper : Fin n) (kl ku : ℝ)
    (hl : 0 ≤ kl) (hu : 0 ≤ ku) (i : Fin n) :
    0 ≤ stationaryWallDiagonal lower upper kl ku i := by
  unfold stationaryWallDiagonal
  exact add_nonneg (by split <;> linarith) (by split <;> linarith)

theorem stationary_wall_force {n : ℕ} (lower upper : Fin n) (kl ku : ℝ)
    (v : Fin n → ℝ) :
    (∑ i, stationaryWallDiagonal lower upper kl ku i * v i) =
      kl * v lower + ku * v upper := by
  simp [stationaryWallDiagonal, add_mul, ite_mul, Finset.sum_add_distrib]

/-- wall sums all stationary boundary conductances at a layer; for one layer
both lower and upper wall coefficients contribute to this same diagonal. -/
def shearOperator {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (conductance : Fin m → ℝ) (wall v : Fin n → ℝ) (i : Fin n) : ℝ :=
  laplace B conductance v i + wall i * v i

def shearDissipation {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (conductance : Fin m → ℝ) (wall v : Fin n → ℝ) : ℝ :=
  (∑ e, conductance e * (gradient B v e)^2) + ∑ i, wall i * v i ^ 2

theorem shear_work {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (conductance : Fin m → ℝ) (wall v : Fin n → ℝ) :
    (∑ i, v i * shearOperator B conductance wall v i) =
      shearDissipation B conductance wall v := by
  unfold shearOperator shearDissipation
  simp only [mul_add, Finset.sum_add_distrib]
  rw [pressure_energy]
  congr 1
  apply Finset.sum_congr rfl
  intro i _
  ring

theorem shear_dissipation_nonnegative {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (conductance : Fin m → ℝ) (wall v : Fin n → ℝ)
    (hk : ∀ e, 0 ≤ conductance e) (hw : ∀ i, 0 ≤ wall i) :
    0 ≤ shearDissipation B conductance wall v := by
  exact add_nonneg
    (Finset.sum_nonneg (fun e _ => mul_nonneg (hk e) (sq_nonneg _)))
    (Finset.sum_nonneg (fun i _ => mul_nonneg (hw i) (sq_nonneg _)))

/-- Exact (M+dt*K)v=M*u+dt*q*V coordinate equations imply the forced work ledger. -/
theorem forced_shear_energy {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (conductance : Fin m → ℝ) (wall volume u v : Fin n → ℝ) (rho dt q : ℝ)
    (step : ∀ i, (rho * volume i) * v i + dt * shearOperator B conductance wall v i =
      (rho * volume i) * u i + dt * (q * volume i)) :
    kineticEnergy (fun i => rho * volume i) v - kineticEnergy (fun i => rho * volume i) u +
      kineticEnergy (fun i => rho * volume i) (fun i => v i - u i) +
      dt * shearDissipation B conductance wall v =
      dt * ∑ i, q * volume i * v i := by
  rw [kinetic_increment, ← shear_work B conductance wall v, Finset.mul_sum,
    ← Finset.sum_add_distrib, Finset.mul_sum]
  apply Finset.sum_congr rfl
  intro i _
  have hi := congrArg (fun x : ℝ => x * v i) (step i)
  dsimp only at hi
  ring_nf at hi ⊢
  linarith only [hi]

/-- Internal interface forces cancel by incidence; only body and wall impulses remain. -/
theorem forced_shear_momentum {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0)
    (conductance : Fin m → ℝ) (wall volume u v : Fin n → ℝ) (rho dt q : ℝ)
    (step : ∀ i, (rho * volume i) * v i + dt * shearOperator B conductance wall v i =
      (rho * volume i) * u i + dt * (q * volume i)) :
    (∑ i, (rho * volume i) * (v i - u i)) =
      dt * (∑ i, q * volume i) - dt * ∑ i, wall i * v i := by
  have internal : (∑ i, laplace B conductance v i) = 0 :=
    internal_flux_conservation B balanced (fun e => conductance e * gradient B v e)
  have coordinate : ∀ i, (rho * volume i) * (v i - u i) =
      dt * (q * volume i) - dt * laplace B conductance v i - dt * (wall i * v i) := by
    intro i
    have hi := step i
    unfold shearOperator at hi
    nlinarith only [hi]
  simp only [coordinate, Finset.sum_sub_distrib, ← Finset.mul_sum, internal,
    mul_zero, sub_zero]

/-- Actual oriented layer edges and lower/upper diagonal give the two-wall ledger. -/
theorem two_wall_shear_momentum {n m : ℕ} (tail head : Fin m → Fin n)
    (lower upper : Fin n) (kl ku : ℝ)
    (conductance : Fin m → ℝ) (volume u v : Fin n → ℝ) (rho dt q : ℝ)
    (step : ∀ i, (rho * volume i) * v i + dt *
      shearOperator (edgeIncidence tail head) conductance
        (stationaryWallDiagonal lower upper kl ku) v i =
      (rho * volume i) * u i + dt * (q * volume i)) :
    (∑ i, (rho * volume i) * (v i - u i)) =
      dt * (∑ i, q * volume i) - dt * (kl * v lower + ku * v upper) := by
  have h := forced_shear_momentum (edgeIncidence tail head) (edge_balanced tail head)
    conductance (stationaryWallDiagonal lower upper kl ku) volume u v rho dt q step
  simpa only [stationary_wall_force] using h

theorem unforced_shear_energy_nonincrease {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (conductance : Fin m → ℝ) (wall volume u v : Fin n → ℝ) (rho dt : ℝ)
    (hr : 0 ≤ rho) (hv : ∀ i, 0 ≤ volume i) (ht : 0 ≤ dt)
    (hk : ∀ e, 0 ≤ conductance e) (hw : ∀ i, 0 ≤ wall i)
    (step : ∀ i, (rho * volume i) * v i + dt * shearOperator B conductance wall v i =
      (rho * volume i) * u i) :
    kineticEnergy (fun i => rho * volume i) v ≤ kineticEnergy (fun i => rho * volume i) u := by
  have he := forced_shear_energy B conductance wall volume u v rho dt 0
    (by simpa using step)
  simp only [zero_mul, Finset.sum_const_zero, mul_zero] at he
  have inc : 0 ≤ kineticEnergy (fun i => rho * volume i) (fun i => v i - u i) := by
    exact Finset.sum_nonneg (fun i _ =>
      mul_nonneg (div_nonneg (mul_nonneg hr (hv i)) (by norm_num)) (sq_nonneg _))
  have diss := mul_nonneg ht (shear_dissipation_nonnegative B conductance wall v hk hw)
  linarith

/-- The gradient of the *specified* fully developed mode u=(u(y),0,0). -/
def reducedShearGradient (derivative : ℝ) (i j : Fin 3) : ℝ :=
  if i = 0 ∧ j = 1 then derivative else 0

/-- mu/2 times the squared symmetric-gradient norm, evaluated at this mode.
This finite algebra does not prove that arbitrary flows have this gradient. -/
theorem reduced_symmetric_strain_dissipation (mu derivative : ℝ) :
    mu / 2 * (∑ i : Fin 3, ∑ j : Fin 3,
      (reducedShearGradient derivative i j + reducedShearGradient derivative j i)^2) =
      mu * derivative^2 := by
  have enumeration : (Finset.univ : Finset (Fin 3)) = {0, 1, 2} := by decide
  simp [enumeration, reducedShearGradient]
  ring

end
end Rheon.ObstacleOperators

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Algebra.BigOperators.Ring.Finset
import Mathlib.Algebra.Order.BigOperators.Group.Finset
import Mathlib.Data.Fintype.Fin
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.FieldSimp

/-! Finite exact-real research contracts. These hypotheses are not established
by mesh queries, a pressure solver, floating-point execution, or a material fit. -/
namespace Rheon.Physics
noncomputable section
open scoped BigOperators

theorem force_work_identity (m u a dt : ℝ) :
    m / 2 * (u + dt * a)^2 - m / 2 * u^2 =
      m * dt * a * u + m / 2 * dt^2 * a^2 := by ring

theorem hydrostatic_increment (u g dt dp rho h : ℝ)
    (hr : rho ≠ 0) (hh : h ≠ 0) (balance : dp = -rho * g * h) :
    u - dt * g - dt * (dp / (rho * h)) = u := by
  rw [balance]
  field_simp
  ring

theorem positive_transmissibility (area rho distance : ℝ)
    (ha : 0 < area) (hr : 0 < rho) (hd : 0 < distance) :
    0 < area / (rho * distance) :=
  div_pos ha (mul_pos hr hd)

theorem weighted_pressure_energy {m : ℕ} (area rho distance gradient : Fin m → ℝ)
    (ha : ∀ e, 0 < area e) (hr : ∀ e, 0 < rho e) (hd : ∀ e, 0 < distance e) :
    0 ≤ ∑ e, (area e / (rho e * distance e)) * (gradient e)^2 := by
  exact Finset.sum_nonneg (fun e _ => mul_nonneg
    (le_of_lt (positive_transmissibility _ _ _ (ha e) (hr e) (hd e))) (sq_nonneg _))

theorem internal_amount_balance {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0) (flux : Fin m → ℝ) :
    (∑ i, ∑ e, B i e * flux e) = 0 := by
  rw [Finset.sum_comm]
  simp_rw [← Finset.sum_mul, balanced, zero_mul]
  simp

theorem sealed_flux_compatibility {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0) (u : Fin m → ℝ)
    (Q source : Fin n → ℝ)
    (constraint : ∀ i, (∑ e, B i e * u e) + Q i = source i) :
    (∑ i, Q i) = ∑ i, source i := by
  have h := Finset.sum_congr (s₁ := Finset.univ) (s₂ := Finset.univ) rfl
    (fun i _ => constraint i)
  rw [Finset.sum_add_distrib, internal_amount_balance B balanced u, zero_add] at h
  exact h

theorem volume_source_balance {n m : ℕ} (B : Fin n → Fin m → ℝ)
    (balanced : ∀ e, (∑ i, B i e) = 0)
    (amount source : Fin n → ℝ) (flux : Fin m → ℝ) (dt : ℝ) :
    (∑ i, (amount i - dt * (∑ e, B i e * flux e) + dt * source i)) =
      (∑ i, amount i) + dt * ∑ i, source i := by
  rw [Finset.sum_add_distrib, Finset.sum_sub_distrib,
    ← Finset.mul_sum, ← Finset.mul_sum, internal_amount_balance B balanced flux]
  ring

theorem weighted_energy_identity {n : ℕ} (mass u v : Fin n → ℝ) :
    (∑ i, mass i * u i ^ 2) - (∑ i, mass i * v i ^ 2) =
      (∑ i, mass i * (u i - v i)^2) +
        2 * ∑ i, mass i * (u i - v i) * v i := by
  rw [← Finset.sum_sub_distrib, Finset.mul_sum, ← Finset.sum_add_distrib]
  apply Finset.sum_congr rfl
  intro i _
  ring

theorem implicit_energy_identity {n : ℕ} (mass u v : Fin n → ℝ)
    (dt dissipation : ℝ)
    (work : (∑ i, mass i * (u i - v i) * v i) = dt * dissipation) :
    (∑ i, mass i * u i ^ 2) - (∑ i, mass i * v i ^ 2) =
      (∑ i, mass i * (u i - v i)^2) + 2 * dt * dissipation := by
  rw [weighted_energy_identity, work]
  ring

theorem implicit_energy_nonincrease {n : ℕ} (mass u v : Fin n → ℝ)
    (dt dissipation : ℝ) (hm : ∀ i, 0 ≤ mass i)
    (ht : 0 ≤ dt) (hd : 0 ≤ dissipation)
    (work : (∑ i, mass i * (u i - v i) * v i) = dt * dissipation) :
    (∑ i, mass i * v i ^ 2) ≤ ∑ i, mass i * u i ^ 2 := by
  have hs : 0 ≤ ∑ i, mass i * (u i - v i)^2 :=
    Finset.sum_nonneg (fun i _ => mul_nonneg (hm i) (sq_nonneg _))
  have hp : 0 ≤ dt * dissipation := mul_nonneg ht hd
  have he := implicit_energy_identity mass u v dt dissipation work
  linarith

theorem slip_power_nonpositive (beta relative_speed : ℝ) (hb : 0 ≤ beta) :
    (-beta * relative_speed) * relative_speed ≤ 0 := by
  nlinarith [mul_nonneg hb (sq_nonneg relative_speed)]

theorem young_adhesion_identity (solid_vapor solid_liquid sigma cosine : ℝ)
    (young : solid_vapor - solid_liquid = sigma * cosine) :
    solid_vapor + sigma - solid_liquid = sigma * (1 + cosine) := by
  nlinarith

theorem adhesion_bounds (sigma cosine : ℝ) (hs : 0 ≤ sigma)
    (hl : -1 ≤ cosine) (hu : cosine ≤ 1) :
    0 ≤ sigma * (1 + cosine) ∧ sigma * (1 + cosine) ≤ 2 * sigma := by
  constructor
  · exact mul_nonneg hs (by linarith)
  · nlinarith [mul_nonneg hs (show 0 ≤ 1 - cosine by linarith)]

end
end Rheon.Physics

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Tactic.Ring
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Positivity
import Mathlib.Tactic.LinearCombination
namespace Rheon.PracticalStudy
noncomputable section
-- Unit gap/area, steady unidirectional Newtonian solution with body force G.
theorem wall_force_balance (mu U G : ℝ) :
 (-mu*U-G/2)+(mu*U-G/2)+G=0 := by ring

theorem wall_power_balance (mu U G : ℝ) (h : mu ≠ 0) :
 U*(mu*U-G/2)+G*(U/2+G/(12*mu))=mu*U^2+G^2/(12*mu) := by
 field_simp
 <;> ring

theorem wall_dissipation_nonnegative (mu U G : ℝ) (h : 0 < mu) :
 0 ≤ mu*U^2+G^2/(12*mu) := by
 exact add_nonneg (mul_nonneg (le_of_lt h) (sq_nonneg U))
   (div_nonneg (sq_nonneg G) (by positivity))

-- Fixed two-unknown Dirichlet matrix K = [[2,-1],[-1,2]].
theorem matrix_energy_identity (x y : ℝ) :
 2*x^2-2*x*y+2*y^2=(x-y)^2+x^2+y^2 := by ring

theorem matrix_energy_positive (x y : ℝ) (h : x ≠ 0 ∨ y ≠ 0) :
 0 < 2*x^2-2*x*y+2*y^2 := by
 rcases h with hx | hy
 · nlinarith [sq_nonneg (x-y), sq_nonneg y, sq_pos_of_ne_zero hx]
 · nlinarith [sq_nonneg (x-y), sq_nonneg x, sq_pos_of_ne_zero hy]

theorem jacobi_smooth_error (x y a : ℝ) :
 (((1-a)*x+a*y/2)+((1-a)*y+a*x/2))/2=(1-a/2)*((x+y)/2) := by ring

theorem jacobi_alternating_error (x y a : ℝ) :
 (((1-a)*x+a*y/2)-((1-a)*y+a*x/2))/2=(1-3*a/2)*((x-y)/2) := by ring

theorem jacobi_factors_strict (a : ℝ) (ha : 0 < a) (hb : a < 4/3) :
 |1-a/2| < 1 ∧ |1-3*a/2| < 1 := by
 constructor <;> apply abs_lt.mpr <;> constructor <;> linarith

theorem jacobi_unstable_mode : (1-3*(3/2:ℝ)/2)^2=25/16 := by norm_num

-- Exact-real arithmetic does not establish a binary32 implementation.
theorem exact_reduction_order (A : ℝ) : (A+1)-A=(A-A)+1 := by ring

-- Fixed masses and one sample. The premise is the mass-weighted grid mean.
theorem transfer_momentum (m v w u du a : ℝ) (h : (1+m)*u=v+m*w) :
 ((1-a)*(u+du)+a*(v+du))+m*((1-a)*(u+du)+a*(w+du))=
 v+m*w+(1+m)*du := by
 linear_combination (1-a)*h

theorem transfer_variance_scale (m x y a : ℝ) :
 (a*x)^2+m*(a*y)^2=a^2*(x^2+m*y^2) := by ring

theorem transfer_variance_nonincrease (V a : ℝ) (hV : 0 ≤ V)
 (ha : 0 ≤ a) (hb : a ≤ 1) : a^2*V ≤ V := by
 have hfactor : 0 ≤ 1-a*a := by nlinarith
 nlinarith [mul_nonneg hfactor hV]

-- Supplied-field constitutive probe; no transport/model-validation theorem.
theorem sst_denominator_positive (a w s f : ℝ) (h : 0 < a*w) :
 0 < max (a*w) (s*f) := lt_of_lt_of_le h (le_max_left _ _)

theorem sst_viscosity_nonnegative (rho a k d : ℝ)
 (hr : 0 ≤ rho) (ha : 0 ≤ a) (hk : 0 ≤ k) (hd : 0 ≤ d) :
 0 ≤ rho*a*k/d := div_nonneg (mul_nonneg (mul_nonneg hr ha) hk) hd

theorem sst_production_bounds (P cap : ℝ) (hP : 0 ≤ P) (hc : 0 ≤ cap) :
 0 ≤ min P cap ∧ min P cap ≤ P ∧ min P cap ≤ cap := by
 exact ⟨le_min hP hc, min_le_left _ _, min_le_right _ _⟩

theorem sst_laminar_viscosity (rho a d : ℝ) : rho*a*0/d=0 := by simp

-- Richardson algebra under an explicitly supplied quadratic error model.
theorem quadratic_extrapolation (q c h : ℝ) :
 (q-c*h^2/4)+((q-c*h^2/4)-(q-c*h^2))/3=q := by ring

theorem constant_bias_survives_refinement (q b c h : ℝ) :
 (q+b-c*h^2/4)+((q+b-c*h^2/4)-(q+b-c*h^2))/3=q+b := by ring

-- Two-cell full residual and gauge invariance, with an explicit tolerance.
theorem candidate_full_residual (e : ℝ) :
 1-(1-e)=e ∧ -1+(1-e)=-e := by constructor <;> ring

theorem candidate_gauge_flux (s e : ℝ) : s-(s-1+e)=1-e := by ring

theorem candidate_divergence_bound (e t : ℝ) (h : |e| ≤ t) :
 |-e| ≤ t ∧ |e| ≤ t := by simpa only [abs_neg] using And.intro h h

theorem positive_tolerance_is_not_zero :
 |(1/16:ℝ)| ≤ 1/8 ∧ (1/16:ℝ) ≠ 0 := by
 constructor
 · rw [abs_of_nonneg (by norm_num : (0:ℝ) ≤ 1/16)]
   norm_num
 · norm_num
end
end Rheon.PracticalStudy

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Tactic.Ring
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.NormNum

/-! Exact two-cell component contracts. Positive physical weights/volumes are
specified in the lesson. These identities do not prove an arbitrary mesh's
nullspace, a PDE solution, floating-point refinement, or a solver's runtime. -/
namespace Rheon.PressureGauge
noncomputable section

theorem component_sum_zero (w p q : ℝ) : w*(p-q)+w*(q-p)=0 := by ring

theorem pair_energy (w p q : ℝ) :
    p*(w*(p-q))+q*(w*(q-p))=w*(q-p)^2 := by ring

theorem component_compatibility (w p q l r : ℝ)
    (hl : w*(p-q)=l) (hr : w*(q-p)=r) : l+r=0 := by
  rw [← hl, ← hr]; ring

theorem independent_gauge_jumps (p q s t a b : ℝ) :
    (q+a)-(p+a)=q-p ∧ (t+b)-(s+b)=t-s := by constructor <;> ring

theorem fixed_cell_residual (w l r : ℝ) (hw : w ≠ 0) :
    l-w*(0-r/w)=l+r ∧ r-w*(r/w-0)=0 := by
  constructor <;> field_simp <;> ring

theorem fixed_cell_solves_compatible (w l r : ℝ) (hw : w ≠ 0)
    (compatible : l+r=0) :
    l-w*(0-r/w)=0 ∧ r-w*(r/w-0)=0 := by
  have h := fixed_cell_residual w l r hw
  exact ⟨h.1.trans compatible, h.2⟩

theorem opposite_component_imbalances :
    (1:ℝ)+0+(-1)+0=0 ∧ (1:ℝ)+0≠0 ∧ (-1:ℝ)+0≠0 := by norm_num

theorem arithmetic_mean_zero (p q : ℝ) :
    (p-(p+q)/2)+(q-(p+q)/2)=0 := by ring

theorem volume_mean_zero (v z p q : ℝ) (hv : v+z ≠ 0) :
    v*(p-(v*p+z*q)/(v+z))+z*(q-(v*p+z*q)/(v+z))=0 := by
  field_simp; ring

theorem projected_component_sum (l r : ℝ) :
    (l-(l+r)/2)+(r-(l+r)/2)=0 := by ring

theorem projected_original_residual (l r : ℝ) :
    l-(l-(l+r)/2)=(l+r)/2 ∧ r-(r-(l+r)/2)=(l+r)/2 := by
  constructor <;> ring

theorem volume_weighted_left_null (v z w p q : ℝ) (hv : v ≠ 0) (hz : z ≠ 0) :
    v*(w*(p-q)/v)+z*(w*(q-p)/z)=0 := by field_simp; ring

theorem boundary_work_shift (p q a qL qR : ℝ) :
    (p+a)*qL+(q+a)*qR-(p*qL+q*qR)=a*(qL+qR) := by ring

theorem balanced_boundary_work (p q a flux : ℝ) :
    (p+a)*flux+(q+a)*(-flux)=p*flux+q*(-flux) := by ring

theorem pressure_pairing_with_boundary (p q flux qL qR : ℝ) :
    p*(flux+qL)+q*(-flux+qR)=-(q-p)*flux+p*qL+q*qR := by ring

end
end Rheon.PressureGauge

Download Lean source

import Rheon.Discrete
import Mathlib.Tactic.Linarith

/-! Exact-real overlap, volume complement and shared-face algebra. These
contracts do not prove Rust mesh admission, IEEE/predicate robustness, component
construction, arbitrary cut meshes or a pressure/viscous solver refinement.
Cells are intervals at represented world endpoints; zero-volume cells are not
pressure unknowns. Active graph weights and component labels are assumptions. -/
namespace Rheon.StaticObstacle
noncomputable section

def overlap (a b c d : ℝ) : ℝ := max 0 (min b d - max a c)

theorem overlap_nonnegative (a b c d : ℝ) : 0 ≤ overlap a b c d := by
  exact le_max_left _ _

theorem overlap_symmetric (a b c d : ℝ) : overlap a b c d = overlap c d a b := by
  simp only [overlap, min_comm b d, max_comm a c]

theorem overlap_bounded (a b c d : ℝ) (interval : a ≤ b) :
    overlap a b c d ≤ b - a := by
  apply max_le
  · linarith
  · have h₁ := min_le_left b d
    have h₂ := le_max_left a c
    linarith

def fluidVolume (hx hy hz dx dy dz : ℝ) : ℝ := hx * hy * hz - dx * dy * dz

theorem fluid_volume_bounds (hx hy hz dx dy dz : ℝ)
    (hx0 : 0 ≤ hx) (hy0 : 0 ≤ hy) (hz0 : 0 ≤ hz)
    (dx0 : 0 ≤ dx) (dy0 : 0 ≤ dy) (dz0 : 0 ≤ dz)
    (dxh : dx ≤ hx) (dyh : dy ≤ hy) (dzh : dz ≤ hz) :
    0 ≤ fluidVolume hx hy hz dx dy dz ∧ fluidVolume hx hy hz dx dy dz ≤ hx * hy * hz := by
  have lower : 0 ≤ dx * dy * dz := mul_nonneg (mul_nonneg dx0 dy0) dz0
  have upper : dx * dy * dz ≤ hx * hy * hz := by
    calc
      dx * dy * dz ≤ hx * dy * dz := mul_le_mul_of_nonneg_right (mul_le_mul_of_nonneg_right dxh dy0) dz0
      _ ≤ hx * hy * dz := mul_le_mul_of_nonneg_right (mul_le_mul_of_nonneg_left dyh hx0) dz0
      _ ≤ hx * hy * hz := mul_le_mul_of_nonneg_left dzh (mul_nonneg hx0 hy0)
  unfold fluidVolume
  constructor <;> linarith

theorem shared_flux_cancels (area speed : ℝ) : area * speed + (-(area * speed)) = 0 := by ring

theorem component_constant_jump_zero {n m : ℕ}
    (tail head : Fin m → Fin n) (weight : Fin m → ℝ)
    (label : Fin n → ℕ) (value : ℕ → ℝ)
    (connected : ∀ e, weight e ≠ 0 → label (tail e) = label (head e)) :
    ∀ e, weight e * (value (label (head e)) - value (label (tail e))) = 0 := by
  intro e
  by_cases zero : weight e = 0
  · simp [zero]
  · rw [connected e zero]
    simp

theorem residual_divergence_scale (volume dt residual divergence : ℝ)
    (active : 0 < volume) (integrated : volume * divergence = dt * residual) :
    divergence = dt * residual / volume := by
  apply (eq_div_iff (ne_of_gt active)).mpr
  simpa only [mul_comm] using integrated

end
end Rheon.StaticObstacle

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Algebra.BigOperators.Ring.Finset
import Mathlib.Algebra.Order.BigOperators.Group.Finset
import Mathlib.Data.Fintype.Fin
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Ring

namespace Rheon.Transport
noncomputable section

def blend (t a b : ℝ) : ℝ := (1 - t) * a + t * b

theorem blend_lower (t a b lo : ℝ) (ht0 : 0 ≤ t) (ht1 : t ≤ 1)
    (ha : lo ≤ a) (hb : lo ≤ b) : lo ≤ blend t a b := by
  unfold blend
  nlinarith [mul_nonneg (sub_nonneg.mpr ht1) (sub_nonneg.mpr ha),
    mul_nonneg ht0 (sub_nonneg.mpr hb)]

theorem blend_upper (t a b hi : ℝ) (ht0 : 0 ≤ t) (ht1 : t ≤ 1)
    (ha : a ≤ hi) (hb : b ≤ hi) : blend t a b ≤ hi := by
  unfold blend
  nlinarith [mul_nonneg (sub_nonneg.mpr ht1) (sub_nonneg.mpr ha),
    mul_nonneg ht0 (sub_nonneg.mpr hb)]

theorem blend_constant (t a : ℝ) : blend t a a = a := by
  unfold blend
  ring

theorem convex_sum_bounds {n : ℕ} (w q : Fin n → ℝ) (lo hi : ℝ)
    (hw : ∀ i, 0 ≤ w i) (hs : (∑ i, w i) = 1)
    (hlo : ∀ i, lo ≤ q i) (hhi : ∀ i, q i ≤ hi) :
    lo ≤ ∑ i, w i * q i ∧ (∑ i, w i * q i) ≤ hi := by
  constructor
  · calc
      lo = ∑ i, w i * lo := by rw [← Finset.sum_mul, hs]; ring
      _ ≤ ∑ i, w i * q i := Finset.sum_le_sum (fun i _ => mul_le_mul_of_nonneg_left (hlo i) (hw i))
  · calc
      (∑ i, w i * q i) ≤ ∑ i, w i * hi := Finset.sum_le_sum (fun i _ => mul_le_mul_of_nonneg_left (hhi i) (hw i))
      _ = hi := by rw [← Finset.sum_mul, hs]; ring

theorem upwind_positive (c qleft q : ℝ) (hc0 : 0 ≤ c) (hc1 : c ≤ 1)
    (hl : 0 ≤ qleft) (hq : 0 ≤ q) : 0 ≤ q - c * (q - qleft) := by
  have h := blend_lower c q qleft 0 hc0 hc1 hq hl
  unfold blend at h
  nlinarith

theorem explicit_diffusion_positive (r qm q qp : ℝ) (hr0 : 0 ≤ r)
    (hr1 : r ≤ 1 / 2) (hm : 0 ≤ qm) (hq : 0 ≤ q) (hp : 0 ≤ qp) :
    0 ≤ q + r * (qm - 2 * q + qp) := by
  nlinarith [mul_nonneg hr0 hm, mul_nonneg hr0 hp,
    mul_nonneg (show 0 ≤ 1 - 2*r by linarith) hq]

theorem cfl_counterexample : (1 : ℚ) - (3/2) * (1 - 0) = -1/2 := by norm_num

theorem interpolation_not_mass_conservative :
    ((1 : ℚ) / 2 * 1 + 1 / 2 * 0) + (1 / 2 * 1 + 1 / 2 * 0) +
      (1 / 2 * 1 + 1 / 2 * 0) ≠ 1 + 0 + 0 := by norm_num

end
end Rheon.Transport

Download Lean source

import Mathlib.Data.Real.Basic
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Ring

/-! Exact unit-cell, one-face liquid-volume contracts. No implementation,
IEEE arithmetic, multidimensional VOF, dynamics or performance theorem. -/
namespace Rheon.VolumeFlux
noncomputable section

/-- Liquid occupies [0,f]; the positive-x swept slab is [1-c,1]. -/
def leftFlux (f c : ℝ) : ℝ := max 0 (f + c - 1)

/-- Liquid occupies [1-f,1]; same fraction, different placement. -/
def rightFlux (f c : ℝ) : ℝ := min f c

theorem left_flux_bounds (f c : ℝ) (hf0 : 0 ≤ f) (hf1 : f ≤ 1)
    (hc0 : 0 ≤ c) (hc1 : c ≤ 1) :
    0 ≤ leftFlux f c ∧ leftFlux f c ≤ f ∧ leftFlux f c ≤ c := by
  unfold leftFlux
  exact ⟨le_max_left _ _, max_le hf0 (by linarith), max_le hc0 (by linarith)⟩

theorem right_flux_bounds (f c : ℝ) (hf0 : 0 ≤ f) (hc0 : 0 ≤ c) :
    0 ≤ rightFlux f c ∧ rightFlux f c ≤ f ∧ rightFlux f c ≤ c := by
  exact ⟨le_min hf0 hc0, min_le_left _ _, min_le_right _ _⟩

/-- One identical oriented amount is subtracted and added; no source/boundary loss. -/
theorem shared_transfer_conserves (donor receiver q : ℝ) :
    (donor - q) + (receiver + q) = donor + receiver := by ring

/-- Availability and receiving capacity are separate obligations. -/
theorem shared_transfer_bounds (donor receiver q : ℝ)
    (hd : donor ≤ 1) (hr : 0 ≤ receiver) (hq : 0 ≤ q)
    (available : q ≤ donor) (capacity : receiver + q ≤ 1) :
    0 ≤ donor - q ∧ donor - q ≤ 1 ∧
    0 ≤ receiver + q ∧ receiver + q ≤ 1 := by
  constructor
  · linarith
  constructor
  · linarith
  constructor
  · linarith
  · exact capacity

theorem same_fraction_different_flux :
    leftFlux (3/10) (1/2) = 0 ∧ rightFlux (3/10) (1/2) = 3/10 := by
  norm_num [leftFlux, rightFlux]

/-- Full-slab transfer from 0.3 to an empty cell, then separately clamp. -/
theorem clamping_creates_volume :
    max (0 : ℝ) (3/10 - 1/2) + min 1 (0 + 1/2) - 3/10 = 1/5 := by
  norm_num

end
end Rheon.VolumeFlux

Download Lean source

import Rheon.Physics

/-! Exact-real contracts for the fixed flat-slab wall-friction update.
These do not certify Rust arithmetic, spatial convergence or wetting. -/
namespace Rheon.WallFriction
noncomputable section
open scoped BigOperators

theorem wall_force_work (beta u wall : ℝ) :
    (-beta * (u - wall)) * u =
      -beta * (u - wall)^2 + (-beta * (u - wall)) * wall := by ring

theorem common_translation (beta u wall shift : ℝ) :
    -beta * ((u + shift) - (wall + shift)) = -beta * (u - wall) := by ring

theorem relative_dissipation_nonnegative {n : ℕ} (beta u wall : Fin n → ℝ)
    (hb : ∀ i, 0 ≤ beta i) : 0 ≤ ∑ i, beta i * (u i - wall i)^2 := by
  exact Finset.sum_nonneg (fun i _ => mul_nonneg (hb i) (sq_nonneg _))

theorem explicit_work_identity {n : ℕ} (mass u delta force : Fin n → ℝ)
    (dt bulk wallDissipation actuatorPower : ℝ)
    (step : ∀ i, mass i * delta i = dt * force i)
    (power : (∑ i, force i * u i) = actuatorPower - bulk - wallDissipation) :
    (∑ i, mass i / 2 * ((u i + delta i)^2 - (u i)^2 - (delta i)^2)) +
      dt * (bulk + wallDissipation - actuatorPower) = 0 := by
  have hp : (∑ i, mass i * delta i * u i) =
      dt * (actuatorPower - bulk - wallDissipation) := by
    calc
      (∑ i, mass i * delta i * u i) = ∑ i, dt * force i * u i := by
        apply Finset.sum_congr rfl
        intro i _
        rw [step]
      _ = dt * ∑ i, force i * u i := by
        rw [Finset.mul_sum]
        apply Finset.sum_congr rfl
        intro i _
        ring
      _ = dt * (actuatorPower - bulk - wallDissipation) := by rw [power]
  calc
    _ = (∑ i, mass i * delta i * u i) +
        dt * (bulk + wallDissipation - actuatorPower) := by
      congr 1
      apply Finset.sum_congr rfl
      intro i _
      ring
    _ = 0 := by rw [hp]; ring

end
end Rheon.WallFriction