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

@@ -15,7 +15,8 @@
         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)
@@ -51,16 +52,13 @@
 
         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

@@ -109,9 +109,9 @@
         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
 
@@ -124,13 +124,13 @@
         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
@@ -142,8 +142,7 @@
         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
@@ -178,7 +177,7 @@
         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)
@@ -238,13 +237,14 @@
         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)
+        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.
+    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)
@@ -253,7 +253,7 @@
         ! 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)
@@ -265,34 +265,30 @@
             ! 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
 
@@ -305,36 +301,49 @@
         ! 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

src/common/include/3dHardcodedIC.fpp

@@ -11,11 +11,11 @@
     real(wp) :: rcut, xcut ! Intermediate variables for creating smooth initial condition
 
     real(wp), dimension(0:n, 0:p) :: rcut_arr
-    integer :: l, q, s ! Iterators for reading input files
-    integer :: start, end ! Ints to keep track of position in file
-    character(len=1000) :: line ! String to store line in ile
-    character(len=25) :: value ! String to store value in line
-    integer :: NJet ! Number of jets
+    integer                       :: l, q, s ! Iterators for reading input files
+    integer                       :: start, end ! Ints to keep track of position in file
+    character(len=1000)           :: line ! String to store line in ile
+    character(len=25)             :: value ! String to store value in line
+    integer                       :: NJet ! Number of jets
 
     eps = 1e-9_wp
 
@@ -31,20 +31,20 @@
 
         open (unit=10, file="jets.csv", status="old", action="read")
         do q = 0, NJet - 1
-            read (10, '(A)') line  ! Read a full line as a string
+            read (10, '(A)') line ! Read a full line as a string
             start = 1
 
             do l = 0, 2
-                end = index(line(start:), ',')  ! Find the next comma
+                end = index(line(start:), ',') ! Find the next comma
                 if (end == 0) then
-                    value = trim(adjustl(line(start:)))  ! Last value in the line
+                    value = trim(adjustl(line(start:))) ! Last value in the line
                 else
-                    value = trim(adjustl(line(start:start + end - 2)))  ! Extract substring
-                    start = start + end  ! Move to next value
+                    value = trim(adjustl(line(start:start + end - 2))) ! Extract substring
+                    start = start + end ! Move to next value
                 end if
                 if (l == 0) then
-                    read (value, *) y_th_arr(q)  ! Convert string to numeric value
-                elseif (l == 1) then
+                    read (value, *) y_th_arr(q) ! Convert string to numeric value
+                else if (l == 1) then
                     read (value, *) z_th_arr(q)
                 else
                     read (value, *) r_th_arr(q)
@@ -184,14 +184,15 @@
         ! This analytic patch used geometry 9
         Mach = 0.1
         if (patch_id == 1) then
-            q_prim_vf(E_idx)%sf(i, j, k) = 101325 + (Mach**2*376.636429464809**2/16)*(cos(2*x_cc(i)/1) + cos(2*y_cc(j)/1))*(cos(2*z_cc(k)/1) + 2)
+            q_prim_vf(E_idx)%sf(i, j, &
+                & k) = 101325 + (Mach**2*376.636429464809**2/16)*(cos(2*x_cc(i)/1) + cos(2*y_cc(j)/1))*(cos(2*z_cc(k)/1) + 2)
             q_prim_vf(momxb + 0)%sf(i, j, k) = Mach*376.636429464809*sin(x_cc(i)/1)*cos(y_cc(j)/1)*sin(z_cc(k)/1)
             q_prim_vf(momxb + 1)%sf(i, j, k) = -Mach*376.636429464809*cos(x_cc(i)/1)*sin(y_cc(j)/1)*sin(z_cc(k)/1)
         end if
 
     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