Integrate ffmt Fortran formatter into MFC

Replace fprettify + indenter.py with ffmt (https://github.com/sbryngelson/ffmt),
a fast, configurable Fortran formatter written in Rust.

Integration:
- toolchain/bootstrap/format.sh: single ffmt call replaces multi-pass
  fprettify + indenter.py loop
- toolchain/pyproject.toml: replace fprettify dependency with ffmt
- ffmt.toml: MFC formatting configuration
- .gitignore: add .ffmt_cache/
- Remove toolchain/indenter.py

Source formatting:
- All .fpp/.f90 files reformatted with ffmt v0.2.0
- Unicode symbols replaced with LaTeX (sigma, pi, partial, etc.)
- Operator spacing, keyword casing, declaration alignment
- Named end statements, comment wrapping at 132 chars
- Doxygen comment block re-wrapping

Install: pip install ffmt
Repo: https://github.com/sbryngelson/ffmt

Author: sbryngelson

.gitignore

@@ -113,3 +113,4 @@ benchmarks/*.png
 cce_*/
 cce_*.log
 run_cce_*.sh
+.ffmt_cache/

ffmt.toml

@@ -0,0 +1,17 @@
+# MFC Fortran formatting configuration
+# These are the defaults — this file makes them explicit.
+
+indent-width = 4
+keyword-case = "lower"
+normalize-keywords = true
+indent-fypp = true
+
+[whitespace]
+relational = true
+logical = true
+plusminus = true
+multdiv = false
+power = false
+assignment = true
+declaration = true
+comma = true

src/common/include/1dHardcodedIC.fpp

@@ -13,23 +13,20 @@
         ! magnetic field
         q_prim_vf(B_idx%end - 1)%sf(i, 0, 0) = 0.1_wp*sin(2._wp*pi*x_cc(i))
         q_prim_vf(B_idx%end)%sf(i, 0, 0) = 0.1_wp*cos(2._wp*pi*x_cc(i))
-
     case (170)
-        ! This hardcoded case can be used to start a simulation with initial conditions given from a known 1D profile (e.g. Cantera, SDtoolbox)
+        ! This hardcoded case can be used to start a simulation with initial conditions given from a known 1D profile (e.g. Cantera,
+        ! SDtoolbox)
         @: HardcodedReadValues()
-
     case (180)
         ! This is patch is hard-coded for test suite optimization used in the
         ! 1D_shuoser cases: "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)"
         if (patch_id == 2) then
             q_prim_vf(contxb + 0)%sf(i, 0, 0) = 1 + 0.2*sin(5*x_cc(i))
         end if
-
     case (181)
         ! This is patch is hard-coded for test suite optimization used in the
         ! 1D_titarevtorro cases: "patch_icpp(2)%alpha_rho(1)": "1 + 0.1*sin(20*x*pi)"
         q_prim_vf(contxb + 0)%sf(i, 0, 0) = 1 + 0.1*sin(20*x_cc(i)*pi)
-
     case (182)
         ! This patch is a hard-coded for test suite optimization (multiple component diffusion)
         x_mid_diffu = 0.05_wp/2.0_wp
@@ -51,16 +48,11 @@
 
         temp = (320.0_wp - 1350.0_wp)*profile_shape + 1350.0_wp
 
-        molar_mass_inv = y1/31.998_wp + &
-                         y2/18.01508_wp + &
-                         y3/16.04256_wp + &
-                         y4/28.0134_wp
+        molar_mass_inv = y1/31.998_wp + y2/18.01508_wp + y3/16.04256_wp + y4/28.0134_wp
 
         q_prim_vf(contxb)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
-
     case default
         call s_int_to_str(patch_id, iStr)
-        call s_mpi_abort("Invalid hcid specified for patch "//trim(iStr))
+        call s_mpi_abort("Invalid hcid specified for patch " // trim(iStr))
     end select
-
 #:enddef

src/common/include/2dHardcodedIC.fpp

@@ -6,7 +6,6 @@
     real(wp) :: factor
     real(wp) :: r0, alpha, r2
     real(wp) :: sinA, cosA
-
     real(wp) :: r_sq
 
     ! # 207
@@ -18,9 +17,7 @@
 #:enddef
 
 #:def Hardcoded2D()
-
     select case (patch_icpp(patch_id)%hcid) ! 2D_hardcoded_ic example case
-
     case (200)
         if (y_cc(j) <= (-x_cc(i)**3 + 1)**(1._wp/3._wp)) then
             ! Volume Fractions
@@ -107,11 +104,10 @@
             pInt = pref + rhoH*9.81_wp*(1.2_wp - intH)
             q_prim_vf(E_idx)%sf(i, j, 0) = pInt + rhoL*9.81_wp*(intH - y_cc(j))
         end if
-
     case (205) ! 2D lung wave interaction problem
-        h = 0.0_wp           !non dim origin y
-        lam = 1.0_wp         !non dim lambda
-        amp = patch_icpp(patch_id)%a(2)         !to be changed later!       !non dim amplitude
+        h = 0.0_wp ! non dim origin y
+        lam = 1.0_wp ! non dim lambda
+        amp = patch_icpp(patch_id)%a(2) ! to be changed later!       !non dim amplitude
 
         intH = amp*sin(2*pi*x_cc(i)/lam - pi/2) + h
 
@@ -122,29 +118,25 @@
             q_prim_vf(advxb)%sf(i, j, 0) = patch_icpp(1)%alpha(1)
             q_prim_vf(advxe)%sf(i, j, 0) = patch_icpp(1)%alpha(2)
         end if
-
     case (206) ! 2D lung wave interaction problem - horizontal domain
-        h = 0.0_wp           !non dim origin y
-        lam = 1.0_wp         !non dim lambda
+        h = 0.0_wp ! non dim origin y
+        lam = 1.0_wp ! non dim lambda
         amp = patch_icpp(patch_id)%a(2)
 
         intL = amp*sin(2*pi*y_cc(j)/lam - pi/2) + h
 
-        if (x_cc(i) > intL) then        !this is the liquid
+        if (x_cc(i) > intL) then ! this is the liquid
             q_prim_vf(contxb)%sf(i, j, 0) = patch_icpp(1)%alpha_rho(1)
             q_prim_vf(contxe)%sf(i, j, 0) = patch_icpp(1)%alpha_rho(2)
             q_prim_vf(E_idx)%sf(i, j, 0) = patch_icpp(1)%pres
             q_prim_vf(advxb)%sf(i, j, 0) = patch_icpp(1)%alpha(1)
             q_prim_vf(advxe)%sf(i, j, 0) = patch_icpp(1)%alpha(2)
         end if
-
     case (207) ! Kelvin Helmholtz Instability
         sigma = 0.05_wp/sqrt(2.0_wp)
         gauss1 = exp(-(y_cc(j) - 0.75_wp)**2/(2.0_wp*sigma**2))
         gauss2 = exp(-(y_cc(j) - 0.25_wp)**2/(2.0_wp*sigma**2))
-        q_prim_vf(momxb + 1)%sf(i, j, 0) = &
-            0.1_wp*sin(4.0_wp*pi*x_cc(i))*(gauss1 + gauss2)
-
+        q_prim_vf(momxb + 1)%sf(i, j, 0) = 0.1_wp*sin(4.0_wp*pi*x_cc(i))*(gauss1 + gauss2)
     case (208) ! Richtmeyer Meshkov Instability
         lam = 1.0_wp
         eps = 1.0e-6_wp
@@ -160,7 +152,6 @@
             q_prim_vf(advxb)%sf(i, j, 0) = alpha_sf6
             q_prim_vf(advxe)%sf(i, j, 0) = alpha_air
         end if
-
     case (250) ! MHD Orszag-Tang vortex
         ! gamma = 5/3
         !   rho = 25/(36*pi)
@@ -173,12 +164,11 @@
 
         q_prim_vf(B_idx%beg)%sf(i, j, 0) = -sin(2._wp*pi*y_cc(j))/sqrt(4._wp*pi)
         q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = sin(4._wp*pi*x_cc(i))/sqrt(4._wp*pi)
-
     case (251) ! RMHD Cylindrical Blast Wave [Mignone, 2006: Section 4.3.1]
         if (x_cc(i)**2 + y_cc(j)**2 < 0.08_wp**2) then
             q_prim_vf(contxb)%sf(i, j, 0) = 0.01
             q_prim_vf(E_idx)%sf(i, j, 0) = 1.0
-        elseif (x_cc(i)**2 + y_cc(j)**2 <= 1._wp**2) then
+        else if (x_cc(i)**2 + y_cc(j)**2 <= 1._wp**2) then
             ! Linear interpolation between r=0.08 and r=1.0
             factor = (1.0_wp - sqrt(x_cc(i)**2 + y_cc(j)**2))/(1.0_wp - 0.08_wp)
             q_prim_vf(contxb)%sf(i, j, 0) = 0.01_wp*factor + 1.e-4_wp*(1.0_wp - factor)
@@ -223,7 +213,6 @@
             q_prim_vf(momxb)%sf(i, j, 0) = -(2._wp/sqrt(r_sq))*(y_cc(j) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
             q_prim_vf(momxb + 1)%sf(i, j, 0) = (2._wp/sqrt(r_sq))*(x_cc(i) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
         end if
-
     case (253) ! MHD Smooth Magnetic Vortex
         ! Section 5.2 of
         ! Implicit hybridized discontinuous Galerkin methods for compressible magnetohydrodynamics
@@ -238,22 +227,21 @@
         q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = x_cc(i)*exp(1 - (x_cc(i)**2 + y_cc(j)**2))/(2.*pi)
 
         ! pressure
-        q_prim_vf(E_idx)%sf(i, j, 0) = 1._wp + (1 - 2._wp*(x_cc(i)**2 + y_cc(j)**2))*exp(1 - (x_cc(i)**2 + y_cc(j)**2))/((2._wp*pi)**3)
-
-    case (260)  ! Gaussian Divergence Pulse
-        !  Bx(x) = 1 + C * erf((x-0.5)/σ)
-        !  ⇒  ∂Bx/∂x = C * (2/√π) * exp[-((x-0.5)/σ)**2] * (1/σ)
-        !  Choose C = ε * σ * √π / 2  ⇒ ∂Bx/∂x = ε * exp[-((x-0.5)/σ)**2]
-        !  ψ is initialized to zero everywhere.
+        q_prim_vf(E_idx)%sf(i, j, &
+                  & 0) = 1._wp + (1 - 2._wp*(x_cc(i)**2 + y_cc(j)**2))*exp(1 - (x_cc(i)**2 + y_cc(j)**2))/((2._wp*pi)**3)
+    case (260) ! Gaussian Divergence Pulse
+        !  Bx(x) = 1 + C * erf((x-0.5)/\sigma)
+        !  =>  \partialBx/\partialx = C * (2/\sqrt\pi) * exp[-((x-0.5)/\sigma)**2] * (1/\sigma)
+        !  Choose C = \epsilon * \sigma * \sqrt\pi / 2  => \partialBx/\partialx = \epsilon * exp[-((x-0.5)/\sigma)**2]
+        !  \psi is initialized to zero everywhere.
 
         eps_mhd = patch_icpp(patch_id)%a(2)
         sigma = patch_icpp(patch_id)%a(3)
         C_mhd = eps_mhd*sigma*sqrt(pi)*0.5_wp
 
         ! B-field
         q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1._wp + C_mhd*erf((x_cc(i) - 0.5_wp)/sigma)
-
-    case (261)  ! Blob
+    case (261) ! Blob
         r0 = 1._wp/sqrt(8._wp)
         r2 = x_cc(i)**2 + y_cc(j)**2
         r = sqrt(r2)
@@ -264,79 +252,80 @@
             ! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp)
             ! q_prim_vf(E_idx)%sf(i,j,0) = 6._wp - q_prim_vf(B_idx%beg)%sf(i,j,0)**2/2._wp
         end if
-
-    case (262)  ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°)
-        ! rotate by α = atan(2)
+    case (262) ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°)
+        ! rotate by \alpha = atan(2)
         alpha = atan(2._wp)
         cosA = cos(alpha)
         sinA = sin(alpha)
         ! projection along shock normal
         r = x_cc(i)*cosA + y_cc(j)*sinA
 
         if (r <= 0.5_wp) then
-            ! LEFT state: ρ=1, v∥=+10, v⊥=0, p=20, B∥=B⊥=5/√(4π)
+            ! LEFT state: \rho=1, v\parallel=+10, v\perp=0, p=20, B\parallel=B\perp=5/\sqrt(4\pi)
             q_prim_vf(contxb)%sf(i, j, 0) = 1._wp
             q_prim_vf(momxb)%sf(i, j, 0) = 10._wp*cosA
             q_prim_vf(momxb + 1)%sf(i, j, 0) = 10._wp*sinA
             q_prim_vf(E_idx)%sf(i, j, 0) = 20._wp
-            q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA &
-                                               - (5._wp/sqrt(4._wp*pi))*sinA
-            q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA &
-                                                   + (5._wp/sqrt(4._wp*pi))*cosA
+            q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA - (5._wp/sqrt(4._wp*pi))*sinA
+            q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA + (5._wp/sqrt(4._wp*pi))*cosA
         else
-            ! RIGHT state: ρ=1, v∥=−10, v⊥=0, p=1, B∥=B⊥=5/√(4π)
+            ! RIGHT state: \rho=1, v\parallel=-10, v\perp=0, p=1, B\parallel=B\perp=5/\sqrt(4\pi)
             q_prim_vf(contxb)%sf(i, j, 0) = 1._wp
             q_prim_vf(momxb)%sf(i, j, 0) = -10._wp*cosA
             q_prim_vf(momxb + 1)%sf(i, j, 0) = -10._wp*sinA
             q_prim_vf(E_idx)%sf(i, j, 0) = 1._wp
-            q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA &
-                                               - (5._wp/sqrt(4._wp*pi))*sinA
-            q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA &
-                                                   + (5._wp/sqrt(4._wp*pi))*cosA
+            q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA - (5._wp/sqrt(4._wp*pi))*sinA
+            q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA + (5._wp/sqrt(4._wp*pi))*cosA
         end if
         ! v^z and B^z remain zero by default
-
     case (270)
         ! This hardcoded case extrudes a 1D profile to initialize a 2D simulation domain
         @: HardcodedReadValues()
-
     case (280)
         ! This is patch is hard-coded for test suite optimization used in the
         ! 2D_isentropicvortex case:
         ! This analytic patch uses geometry 2
         if (patch_id == 1) then
-            q_prim_vf(E_idx)%sf(i, j, 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*pi))*(5.0/(8.0*1.0*(1.4 + 1.0)*pi))*exp(2.0*1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**(1.4 + 1.0)
-            q_prim_vf(contxb + 0)%sf(i, j, 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*pi))*(5.0/(8.0*1.0*(1.4 + 1.0)*pi))*exp(2.0*1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**1.4
-            q_prim_vf(momxb + 0)%sf(i, j, 0) = 0.0 + (y_cc(j) - patch_icpp(1)%y_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0))
-            q_prim_vf(momxb + 1)%sf(i, j, 0) = 0.0 - (x_cc(i) - patch_icpp(1)%x_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0))
+            q_prim_vf(E_idx)%sf(i, j, &
+                      & 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*pi))*(5.0/(8.0*1.0*(1.4 + 1.0)*pi))*exp(2.0*1.0*(1.0 - (x_cc(i) &
+                      & - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**(1.4 + 1.0)
+            q_prim_vf(contxb + 0)%sf(i, j, &
+                      & 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*pi))*(5.0/(8.0*1.0*(1.4 + 1.0)*pi))*exp(2.0*1.0*(1.0 - (x_cc(i) &
+                      & - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**1.4
+            q_prim_vf(momxb + 0)%sf(i, j, &
+                      & 0) = 0.0 + (y_cc(j) - patch_icpp(1)%y_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1) &
+                      & %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0))
+            q_prim_vf(momxb + 1)%sf(i, j, &
+                      & 0) = 0.0 - (x_cc(i) - patch_icpp(1)%x_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1) &
+                      & %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0))
         end if
-
     case (281)
         ! This is patch is hard-coded for test suite optimization used in the
         ! 2D_acoustic_pulse case:
         ! This analytic patch uses geometry 2
         if (patch_id == 2) then
-            q_prim_vf(E_idx)%sf(i, j, 0) = 101325*(1 - 0.5*(1.4 - 1)*(0.4)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1.4/(1.4 - 1))
-            q_prim_vf(contxb + 0)%sf(i, j, 0) = 1*(1 - 0.5*(1.4 - 1)*(0.4)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1/(1.4 - 1))
+            q_prim_vf(E_idx)%sf(i, j, &
+                      & 0) = 101325*(1 - 0.5*(1.4 - 1)*(0.4)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1.4/(1.4 - 1))
+            q_prim_vf(contxb + 0)%sf(i, j, &
+                      & 0) = 1*(1 - 0.5*(1.4 - 1)*(0.4)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1/(1.4 - 1))
         end if
-
     case (282)
         ! This is patch is hard-coded for test suite optimization used in the
         ! 2D_zero_circ_vortex case:
         ! This analytic patch uses geometry 2
         if (patch_id == 2) then
-            q_prim_vf(E_idx)%sf(i, j, 0) = 101325*(1 - 0.5*(1.4 - 1)*(0.1/0.3)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1.4/(1.4 - 1))
-            q_prim_vf(contxb + 0)%sf(i, j, 0) = 1*(1 - 0.5*(1.4 - 1)*(0.1/0.3)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1/(1.4 - 1))
-            q_prim_vf(momxb + 0)%sf(i, j, 0) = 112.99092883944267*(1 - (0.1/0.3))*y_cc(j)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2)))
+            q_prim_vf(E_idx)%sf(i, j, &
+                      & 0) = 101325*(1 - 0.5*(1.4 - 1)*(0.1/0.3)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1.4/(1.4 - 1))
+            q_prim_vf(contxb + 0)%sf(i, j, &
+                      & 0) = 1*(1 - 0.5*(1.4 - 1)*(0.1/0.3)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1/(1.4 - 1))
+            q_prim_vf(momxb + 0)%sf(i, j, &
+                      & 0) = 112.99092883944267*(1 - (0.1/0.3))*y_cc(j)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2)))
             q_prim_vf(momxb + 1)%sf(i, j, 0) = 112.99092883944267*((0.1/0.3))*x_cc(i)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2)))
         end if
-
     case default
         if (proc_rank == 0) then
             call s_int_to_str(patch_id, iStr)
-            call s_mpi_abort("Invalid hcid specified for patch "//trim(iStr))
+            call s_mpi_abort("Invalid hcid specified for patch " // trim(iStr))
         end if
-
     end select
-
 #:enddef