From 0fbedeef8eb18f46bb022c31bd76a03f6e2a75dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 15:24:50 -0600 Subject: [PATCH 01/18] fix: scalar divergence assignment (cherry picked from commit 44976dcbf011cb2564263bafd2f4bccbffc8ba64) --- src/formal/scalar_2D_s.F90 | 8 +- test/divergence_assignment_2D_test_m.F90 | 167 +++++++++++++++++++++++ test/driver.f90 | 2 + 3 files changed, 174 insertions(+), 3 deletions(-) create mode 100644 test/divergence_assignment_2D_test_m.F90 diff --git a/src/formal/scalar_2D_s.F90 b/src/formal/scalar_2D_s.F90 index c419e66..51215c0 100644 --- a/src/formal/scalar_2D_s.F90 +++ b/src/formal/scalar_2D_s.F90 @@ -15,6 +15,7 @@ ,operator(.isAtLeast.) & ,string_t use tensors_1D_m, only : cell_centers_extended_1D, scalar_1D_t + use differential_operators_1D_m, only : gradient_operator_1D_t use interpolator_1D_m, only : centers_to_faces_1D_t implicit none @@ -195,10 +196,10 @@ allocate(lhs%points_(1,1,1,1)%values_(rhs%cells_(x_dir)+2, rhs%cells_(y_dir)+2)) associate( & - x_last => size(rhs%points_(1,1,1,1)%values_,x_dir) - 1 & - ,y_last => size(rhs%points_(1,1,1,1)%values_,y_dir) - 1 & + x_last => size(rhs%points_(1,1,1,1)%values_,x_dir) + 2 & + ,y_last => size(rhs%points_(1,1,1,1)%values_,y_dir) + 2 & ) - lhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) = rhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) ! internal points + lhs%points_(1,1,1,1)%values_(2:x_last-1, 2:y_last-1) = rhs%points_(1,1,1,1)%values_ ! internal points lhs%points_(1,1,1,1)%values_(1 , : ) = 0E0 ! x_min boundary lhs%points_(1,1,1,1)%values_( x_last , : ) = 0E0 ! x_max boundary lhs%points_(1,1,1,1)%values_( : , 1 ) = 0E0 ! y_min boundary @@ -209,6 +210,7 @@ lhs%x_min_ = rhs%x_min_ lhs%x_max_ = rhs%x_max_ lhs%order_ = rhs%order_ + lhs%gradient_operator_1D_ = gradient_operator_1D_t(k=lhs%order_, dx=(lhs%x_max_ - lhs%x_min_)/lhs%cells_, cells=lhs%cells_) call_julienne_assert(lhs%consistent()) call_julienne_assert(lhs%conformable(rhs)) diff --git a/test/divergence_assignment_2D_test_m.F90 b/test/divergence_assignment_2D_test_m.F90 new file mode 100644 index 0000000..eae2b49 --- /dev/null +++ b/test/divergence_assignment_2D_test_m.F90 @@ -0,0 +1,167 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module divergence_assignment_2D_test_m + use julienne_m, only : & + operator(//) & + ,operator(.all.) & + ,operator(.also.) & + ,operator(.equalsExpected.) & + ,operator(.approximates.) & + ,operator(.within.) & + ,passing_test & + ,string_t & + ,test_description_t & + ,test_diagnosis_t & + ,test_result_t & + ,test_t & + ,usher + use formal_m, only : scalar_2D_t, vector_2D_t, gradient_2D_t, divergence_2D_t, scalar_2D_initializer_i + implicit none + + integer, parameter :: x_cells_1 = 9, y_cells_1 = 10 + integer, parameter :: x_cells_2 = 11, y_cells_2 = 12 + real, parameter :: tolerance = 1E-6 + + type, extends(test_t) :: divergence_assignment_2D_test_t + contains + procedure, nopass :: subject + procedure, nopass :: results + end type + +contains + + pure function subject() result(test_subject) + character(len=:), allocatable :: test_subject + test_subject = 'The scalar_2D_t divergence assignment' + end function + + function results() result(test_results) + type(divergence_assignment_2D_test_t) test_case + type(test_result_t), allocatable :: test_results(:) + + test_results = test_case%run([ & + test_description_t('copying a nonsquare divergence field into a scalar, preserving boundaries, repeated assignment, and .grad.', usher(check_divergence_assignment)) & + ]) + end function + + pure function pattern_a(x, y) result(values) + real, intent(in) :: x(:), y(:) + real values(size(x), size(y)) + do concurrent(integer :: i = 1:size(x), j = 1:size(y)) default(none) shared(x, y, values) + values(i,j) = 1000E0 + 10E0*x(i) + y(j) + end do + end function + + pure function pattern_b(x, y) result(values) + real, intent(in) :: x(:), y(:) + real values(size(x), size(y)) + do concurrent(integer :: i = 1:size(x), j = 1:size(y)) default(none) shared(x, y, values) + values(i,j) = -500E0 + x(i) - 3E0*y(j) + end do + end function + + function check_divergence_assignment() result(test_diagnosis) + type(test_diagnosis_t) test_diagnosis + test_diagnosis = passing_test() + test_diagnosis = test_diagnosis .also. check_case(4, [x_cells_1, y_cells_1], [-2E0,-1E0], [1E0,3E0]) + test_diagnosis = test_diagnosis .also. check_case(2, [x_cells_2, y_cells_2], [-1E0,-2E0], [2E0,4E0]) + end function + + function check_case(order, cells, x_min, x_max) result(test_diagnosis) + integer, intent(in) :: order, cells(2) + real, intent(in) :: x_min(2), x_max(2) + type(test_diagnosis_t) test_diagnosis + procedure(scalar_2D_initializer_i), pointer :: initializer_a, initializer_b + type(divergence_2D_t) rhs_a, rhs_b + type(scalar_2D_t) lhs, expected_a, expected_b + type(gradient_2D_t) grad_lhs, grad_expected + real, allocatable :: lhs_values(:,:), rhs_values(:,:) + real, allocatable :: expected_values(:,:) + + test_diagnosis = passing_test() + initializer_a => pattern_a + initializer_b => pattern_b + + rhs_a = divergence_2D_t(initializer_a, order=order, cells=cells, x_min=x_min, x_max=x_max) + rhs_b = divergence_2D_t(initializer_b, order=order, cells=cells, x_min=x_min, x_max=x_max) + allocate(expected_values(cells(1)+2, cells(2)+2)) + expected_values = 0E0 + expected_values(2:size(expected_values,1)-1, 2:size(expected_values,2)-1) = rhs_a%values() + expected_a = scalar_2D_t(expected_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) + expected_values = 0E0 + expected_values(2:size(expected_values,1)-1, 2:size(expected_values,2)-1) = rhs_b%values() + expected_b = scalar_2D_t(expected_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) + + lhs = rhs_a + lhs_values = lhs%values() + rhs_values = rhs_a%values() + test_diagnosis = test_diagnosis .also. (check_shape(lhs_values)) .also. (check_contents(lhs_values, rhs_values)) .also. (check_gradient(lhs, expected_a)) + + lhs = rhs_b + lhs_values = lhs%values() + rhs_values = rhs_b%values() + test_diagnosis = test_diagnosis .also. (check_shape(lhs_values)) .also. (check_contents(lhs_values, rhs_values)) .also. (check_gradient(lhs, expected_b)) + + lhs = rhs_a + lhs_values = lhs%values() + rhs_values = rhs_a%values() + test_diagnosis = test_diagnosis .also. (check_shape(lhs_values)) .also. (check_contents(lhs_values, rhs_values)) .also. (check_gradient(lhs, expected_a)) + + contains + + pure function expected_initializer(x, y) result(values) + real, intent(in) :: x(:), y(:) + real values(size(x), size(y)) + values = expected_values + end function + + function check_shape(values) result(diagnosis) + real, intent(in) :: values(:,:) + type(test_diagnosis_t) diagnosis + diagnosis = passing_test() + diagnosis = diagnosis .also. (.all. (shape(values) .equalsExpected. [cells(1) + 2, cells(2) + 2])) + end function + + function check_contents(values, rhs_values) result(diagnosis) + real, intent(in) :: values(:,:), rhs_values(:,:) + type(test_diagnosis_t) diagnosis + diagnosis = passing_test() + diagnosis = diagnosis .also. (boundary_zero(values)) .also. (interior_copy(values, rhs_values)) + end function + + function boundary_zero(values) result(diagnosis) + real, intent(in) :: values(:,:) + type(test_diagnosis_t) diagnosis + integer :: nx, ny + nx = size(values,1) + ny = size(values,2) + diagnosis = passing_test() + diagnosis = diagnosis .also. (maxval(abs(values(1,:))) .approximates. 0E0 .within. tolerance) + diagnosis = diagnosis .also. (maxval(abs(values(nx,:))) .approximates. 0E0 .within. tolerance) + diagnosis = diagnosis .also. (maxval(abs(values(:,1))) .approximates. 0E0 .within. tolerance) + diagnosis = diagnosis .also. (maxval(abs(values(:,ny))) .approximates. 0E0 .within. tolerance) + end function + + function interior_copy(values, rhs_values) result(diagnosis) + real, intent(in) :: values(:,:), rhs_values(:,:) + type(test_diagnosis_t) diagnosis + diagnosis = passing_test() + diagnosis = diagnosis .also. (.all. (values(2:size(values,1)-1, 2:size(values,2)-1) .approximates. rhs_values .within. tolerance)) + end function + + function check_gradient(scalar_field, expected_scalar) result(diagnosis) + type(scalar_2D_t), intent(in) :: scalar_field + type(scalar_2D_t), intent(in) :: expected_scalar + type(test_diagnosis_t) diagnosis + type(gradient_2D_t) grad_expected + diagnosis = passing_test() + grad_lhs = .grad. scalar_field + grad_expected = .grad. expected_scalar + diagnosis = diagnosis .also. (.all. (shape(grad_lhs%to_centers_extended()) .equalsExpected. shape(grad_expected%to_centers_extended()))) + diagnosis = diagnosis .also. (.all. (grad_lhs%to_centers_extended() .approximates. grad_expected%to_centers_extended() .within. tolerance)) + end function + + end function check_case + +end module divergence_assignment_2D_test_m diff --git a/test/driver.f90 b/test/driver.f90 index 35d7202..2f7dd95 100644 --- a/test/driver.f90 +++ b/test/driver.f90 @@ -11,6 +11,7 @@ program test_suite_driver use interpolator_1D_test_m, only : interpolator_1D_test_t use scalar_1D_test_m, only : scalar_1D_test_t use scalar_2D_test_m, only : scalar_2D_test_t + use divergence_assignment_2D_test_m, only : divergence_assignment_2D_test_t use scalar_3D_test_m, only : scalar_3D_test_t use vector_2D_test_m, only : vector_2D_test_t use vector_3D_test_m, only : vector_3D_test_t @@ -24,6 +25,7 @@ program test_suite_driver ,test_fixture_t(interpolator_1D_test_t()) & ,test_fixture_t(scalar_1D_test_t()) & ,test_fixture_t(scalar_2D_test_t()) & + ,test_fixture_t(divergence_assignment_2D_test_t()) & ,test_fixture_t(scalar_3D_test_t()) & ,test_fixture_t(vector_2D_test_t()) & ,test_fixture_t(vector_3D_test_t()) & From 363f8999e7c4575d8926c669fbb2ce8eab7cea06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 21:25:30 -0600 Subject: [PATCH 02/18] test: use portable scalar initializer callbacks Pass explicit procedure pointers to the existing scalar constructors. Use pure module procedures for expected fields so Flang on macOS does not require executable-stack trampolines for a host-associated callback. Preserve the zero-boundary oracle and all numerical assertions. --- test/divergence_assignment_2D_test_m.F90 | 33 ++++++++++++++---------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/test/divergence_assignment_2D_test_m.F90 b/test/divergence_assignment_2D_test_m.F90 index eae2b49..83cb342 100644 --- a/test/divergence_assignment_2D_test_m.F90 +++ b/test/divergence_assignment_2D_test_m.F90 @@ -61,6 +61,20 @@ pure function pattern_b(x, y) result(values) end do end function + pure function expected_pattern_a(x, y) result(values) + real, intent(in) :: x(:), y(:) + real values(size(x), size(y)) + values = 0E0 + values(2:size(x)-1, 2:size(y)-1) = pattern_a(x(2:size(x)-1), y(2:size(y)-1)) + end function + + pure function expected_pattern_b(x, y) result(values) + real, intent(in) :: x(:), y(:) + real values(size(x), size(y)) + values = 0E0 + values(2:size(x)-1, 2:size(y)-1) = pattern_b(x(2:size(x)-1), y(2:size(y)-1)) + end function + function check_divergence_assignment() result(test_diagnosis) type(test_diagnosis_t) test_diagnosis test_diagnosis = passing_test() @@ -73,25 +87,22 @@ function check_case(order, cells, x_min, x_max) result(test_diagnosis) real, intent(in) :: x_min(2), x_max(2) type(test_diagnosis_t) test_diagnosis procedure(scalar_2D_initializer_i), pointer :: initializer_a, initializer_b + procedure(scalar_2D_initializer_i), pointer :: expected_initializer_a, expected_initializer_b type(divergence_2D_t) rhs_a, rhs_b type(scalar_2D_t) lhs, expected_a, expected_b type(gradient_2D_t) grad_lhs, grad_expected real, allocatable :: lhs_values(:,:), rhs_values(:,:) - real, allocatable :: expected_values(:,:) test_diagnosis = passing_test() initializer_a => pattern_a initializer_b => pattern_b + expected_initializer_a => expected_pattern_a + expected_initializer_b => expected_pattern_b rhs_a = divergence_2D_t(initializer_a, order=order, cells=cells, x_min=x_min, x_max=x_max) rhs_b = divergence_2D_t(initializer_b, order=order, cells=cells, x_min=x_min, x_max=x_max) - allocate(expected_values(cells(1)+2, cells(2)+2)) - expected_values = 0E0 - expected_values(2:size(expected_values,1)-1, 2:size(expected_values,2)-1) = rhs_a%values() - expected_a = scalar_2D_t(expected_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) - expected_values = 0E0 - expected_values(2:size(expected_values,1)-1, 2:size(expected_values,2)-1) = rhs_b%values() - expected_b = scalar_2D_t(expected_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) + expected_a = scalar_2D_t(expected_initializer_a, order=order, cells=cells, x_min=x_min, x_max=x_max) + expected_b = scalar_2D_t(expected_initializer_b, order=order, cells=cells, x_min=x_min, x_max=x_max) lhs = rhs_a lhs_values = lhs%values() @@ -110,12 +121,6 @@ function check_case(order, cells, x_min, x_max) result(test_diagnosis) contains - pure function expected_initializer(x, y) result(values) - real, intent(in) :: x(:), y(:) - real values(size(x), size(y)) - values = expected_values - end function - function check_shape(values) result(diagnosis) real, intent(in) :: values(:,:) type(test_diagnosis_t) diagnosis From fdd622f03eea8e9ef2bf173626f4bbd6990e1213 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 15:24:53 -0600 Subject: [PATCH 03/18] perf: fold flip permutation into one assignment (cherry picked from commit 38661318d2afaaf5290bf1b0a9d2734154ec45f1) --- src/formal/differential_operators_1D_m.F90 | 11 +---------- src/formal/divergence_operator_1D_s.F90 | 12 +----------- src/formal/gradient_operator_1D_s.F90 | 12 +----------- 3 files changed, 3 insertions(+), 32 deletions(-) diff --git a/src/formal/differential_operators_1D_m.F90 b/src/formal/differential_operators_1D_m.F90 index b79f3bf..36c764f 100644 --- a/src/formal/differential_operators_1D_m.F90 +++ b/src/formal/differential_operators_1D_m.F90 @@ -145,16 +145,7 @@ pure function negate_and_flip(A) result(Ap) real, allocatable :: Ap(:,:) allocate(Ap, mold=A) - - reverse_elements_within_rows_and_flip_sign: & - do concurrent(integer :: row = 1:size(Ap,1)) default(none) shared(Ap, A) - Ap(row,:) = -A(row,size(A,2):1:-1) - end do reverse_elements_within_rows_and_flip_sign - - reverse_elements_within_columns: & - do concurrent(integer :: column = 1 : size(Ap,2)) default(none) shared(Ap) - Ap(:,column) = Ap(size(Ap,1):1:-1,column) - end do reverse_elements_within_columns + Ap = -A(size(A,1):1:-1,size(A,2):1:-1) end function diff --git a/src/formal/divergence_operator_1D_s.F90 b/src/formal/divergence_operator_1D_s.F90 index c9cf6dd..a4170e2 100644 --- a/src/formal/divergence_operator_1D_s.F90 +++ b/src/formal/divergence_operator_1D_s.F90 @@ -18,19 +18,9 @@ pure function negate_and_flip(A) result(Ap) !! Transform a mimetic matrix upper block into a lower block real, intent(in) :: A(:,:) real, allocatable :: Ap(:,:) - integer row, column allocate(Ap, mold=A) - - reverse_elements_within_rows_and_flip_sign: & - do concurrent(row = 1:size(Ap,1)) - Ap(row,:) = -A(row,size(A,2):1:-1) - end do reverse_elements_within_rows_and_flip_sign - - reverse_elements_within_columns: & - do concurrent(column = 1 : size(Ap,2)) - Ap(:,column) = Ap(size(Ap,1):1:-1,column) - end do reverse_elements_within_columns + Ap = -A(size(A,1):1:-1,size(A,2):1:-1) end function diff --git a/src/formal/gradient_operator_1D_s.F90 b/src/formal/gradient_operator_1D_s.F90 index 673f50b..435b56c 100644 --- a/src/formal/gradient_operator_1D_s.F90 +++ b/src/formal/gradient_operator_1D_s.F90 @@ -19,19 +19,9 @@ pure function negate_and_flip(A) result(Ap) !! Transform a mimetic matrix upper block into a lower block real, intent(in) :: A(:,:) real, allocatable :: Ap(:,:) - integer row, column allocate(Ap, mold=A) - - reverse_elements_within_rows_and_flip_sign: & - do concurrent(row = 1:size(Ap,1)) - Ap(row,:) = -A(row,size(A,2):1:-1) - end do reverse_elements_within_rows_and_flip_sign - - reverse_elements_within_columns: & - do concurrent(column = 1 : size(Ap,2)) - Ap(:,column) = Ap(size(Ap,1):1:-1,column) - end do reverse_elements_within_columns + Ap = -A(size(A,1):1:-1,size(A,2):1:-1) end function From 4d3d91789723b6de45e96326dd105c517da0092c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 18:25:32 -0600 Subject: [PATCH 04/18] perf: keep advection RK4 workspace resident on the GPU Preserve the operator-based solver as a numerical reference, retain normal CSV output, and add substantial benchmark and raw-field modes. On M4 Max with the coalescing LFortran build, the 2048x2048 Gaussian case takes 2.196 s wall / 2.186 s GPU and matches native Metal within 5.96e-8. (cherry picked from commit 39defc6ae58d028e6321359868d5e90e5dea2090) --- README.md | 24 + example/2D-advection-diffusion.F90 | 480 +++++++++++- example/formal_advection_workspace.F90 | 684 ++++++++++++++++++ example/gpu_benchmark_m.F90 | 123 ++++ example/gpu_benchmark_shim.c | 158 ++++ example/scripts/compare_2d_fields.py | 221 ++++++ .../validate_2d_advection_workspace.py | 142 ++++ 7 files changed, 1798 insertions(+), 34 deletions(-) create mode 100644 example/formal_advection_workspace.F90 create mode 100644 example/gpu_benchmark_m.F90 create mode 100644 example/gpu_benchmark_shim.c create mode 100755 example/scripts/compare_2d_fields.py create mode 100644 example/scripts/validate_2d_advection_workspace.py diff --git a/README.md b/README.md index 04e1b14..e748482 100644 --- a/README.md +++ b/README.md @@ -69,6 +69,30 @@ where `s` is the concentration of a passive scalar quantity, `D` is a molecular coefficient, and `v` is a prescribed velocity field. This repository's 2D-advection-diffusion [program](./example/2D-advection-diffusion.F90), for example, demonstrates how to advance the above equation in time using a Runge-Kutta scheme. +It accepts `--solver workspace|operators`, `--cells`, `--steps`, `--dt`, and +`--order` flags for benchmarking; without arguments it retains the 51×51, +500-step, `dt=1e-4` regression case, and when `--cells` changes without `--dt` +the default time step is scaled by the square of the grid refinement factor. +Use `--benchmark` for the 257×257, 30000-step preset, `--no-csv` to suppress +CSV output on large runs, `--dump-input ` to write the native FADV2D1 raw +input snapshot before warmup/timing, and `--dump-final-raw ` for binary +final-field output. The tracked comparison helper lives at +[`example/scripts/compare_2d_fields.py`](./example/scripts/compare_2d_fields.py). +For GPU benchmarking, compile with `-DFORMAL_ENABLE_GPU_BENCHMARK`, pass +`--gpu-profile`, and set `LFORTRAN_GPU_PROFILE=1` to collect Metal telemetry +through the runtime snapshot shim. + +Build LFortran with Pixi, then run its application benchmark entry point. +It supplies the C-shim flags, checks the actual Metal executable, and runs +it under the GPU guard: +```bash +export FORMAL_ROOT=/path/to/formal +cd /path/to/lfortran +pixi run -e llvm11 bash examples/formal_advection/run_formal.sh +``` +The default is the measured 2048 x 2048, 1000-step case; append example +arguments to select another workload. Build directories are keyed by the +compiler binary hash to avoid stale kernels after compiler changes. ### Other example programs See this repository's [example](./example) subdirectory for additional demonstrations of diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90 index 69a474c..79cef80 100644 --- a/example/2D-advection-diffusion.F90 +++ b/example/2D-advection-diffusion.F90 @@ -33,54 +33,262 @@ pure function taylor_green_velocity(x,y) result(velocity) program advection_diffusion_2D !! Solve the advection-diffusion equation for a passive scalar moving through a !! static velocity field define by 2D Taylor-Green vortices. + use formal_advection_workspace, only : & + fa_advance_rk4 & + ,fa_fields_t & + ,fa_init_from_arrays & + ,fa_read_raw_input & + ,fa_validate_cfg & + ,fa_validate_working_set & + ,fa_reconstruct_scalar & + ,fa_write_final_raw use fields_m, only : scalar_field, taylor_green_velocity - use formal_m, only : scalar_2D_t, vector_2D_t, scalar_2D_initializer_i, vector_2D_initializer_i + use gpu_benchmark_m, only : gpu_benchmark_begin, gpu_benchmark_end, gpu_benchmark_stats_t + use formal_m, only : scalar_2D_t, vector_2D_t, scalar_2D_initializer_i, vector_2D_initializer_i, x_dir, y_dir + use iso_fortran_env, only : int64, real32, real64, output_unit + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite + use julienne_m, only : command_line_t implicit none - integer, parameter :: order = 4 + integer, parameter :: default_order = 4 + integer, parameter :: solver_workspace = 1, solver_operators = 2 + real, parameter :: default_dt = 1E-4 + real, parameter :: pi = acos(-1E0) + real, parameter :: x_min(2) = [-pi,-pi], x_max(2) = [pi,pi] + type numerical_arguments_t + integer :: order = default_order + integer :: steps = 500 + real :: dt = default_dt + integer :: cells(2) = [51,51] + integer :: solver = 0 + logical :: benchmark = .false. + logical :: no_csv = .false. + logical :: dt_supplied = .false. + logical :: cells_supplied = .false. + logical :: steps_supplied = .false. + character(len=:), allocatable :: input_path + character(len=:), allocatable :: dump_final_raw_path + end type + type(command_line_t) command_line + type(numerical_arguments_t) args + type(gpu_benchmark_stats_t) warmup_gpu_stats, solve_gpu_stats + type(fa_fields_t) workspace_fields procedure(scalar_2D_initializer_i), pointer :: scalar_2D_initializer procedure(vector_2D_initializer_i), pointer :: velocity_2D_initializer - type(scalar_2D_t) s - real, parameter :: pi = acos(-1E0) + type(scalar_2D_t) s, s_initial, s_final, warmup_rhs + type(vector_2D_t) v + real(real32), allocatable :: workspace_initial(:,:) + integer(int64) :: clock_rate, setup_start, setup_end, initial_io_start, initial_io_end + integer(int64) :: solve_start, solve_end, final_io_start, final_io_end + real(real64) :: setup_seconds, initial_io_seconds, solve_seconds, final_io_seconds + real(real64) :: initial_checksum, final_checksum, l_inf_change, relative_l2_change, warmup_checksum, warmup_jit_seconds + real(real64) :: initial_norm_squared, difference_norm_squared + logical :: initial_all_finite, final_all_finite + logical :: gpu_profile_requested + character(len=:), allocatable :: dump_input_path + integer :: step + + if (command_line%argument_present([character(len=len("--help")) :: ("--help"), "-h"])) then + stop new_line('') // new_line('') & + // 'Usage:' // new_line('') & + // ' fpm run \' // new_line('') & + // ' --example 2D-advection-diffusion \' // new_line('') & + // ' --compiler flang \' // new_line('') & + // ' --flag "-O3" \' // new_line('') & + // ' -- [--solver workspace|operators] [--cells ] [--steps ] [--dt ] [--order ] [--benchmark] [--no-csv] [--dump-input ] [--dump-final-raw ] [--input ] [--gpu-profile]' // new_line('') & + // ' or use the no-argument legacy regression case (51x51, 500 steps, dt = 1e-4).' // new_line('') // new_line('') & + // 'When --dt is omitted, the default time step scales as 1e-4 * min(51/Nx, 51/Ny)^2.' // new_line('') // new_line('') & + // 'Use --benchmark for the 257x257, 30000-step preset; explicit --cells, --steps, and --dt override the preset.' // new_line('') // new_line('') & + // 'The workspace solver is the default for order 4; order 2 defaults to operators and rejects workspace explicitly.' // new_line('') // new_line('') & + // 'When --dump-input is set, the initial state is written before warmup/timing in the native FADV2D1 raw format.' // new_line('') // new_line('') & + // 'When --dump-final-raw is set, the final scalar field is written in the same raw format for binary comparison.' // new_line('') // new_line('') & + // 'When --input is set, raw input mode is intended for workspace validation with --no-csv.' // new_line('') // new_line('') & + // 'When --gpu-profile is enabled, set LFORTRAN_GPU_PROFILE=1 to collect Metal telemetry.' // new_line('') // new_line('') & + // 'Examples:' // new_line('') & + // ' fpm run --example 2D-advection-diffusion -- --cells 129,129 --steps 2000 --dt 2e-5' // new_line('') & + // ' fpm run --example 2D-advection-diffusion -- --cells 51,51 --steps 500 --dt 1e-4 --gpu-profile' // new_line('') + end if + + args = get_numerical_arguments() + dump_input_path = command_line%flag_value("--dump-input") + if (.not. allocated(dump_input_path)) dump_input_path = "" + if (.not. allocated(args%dump_final_raw_path)) args%dump_final_raw_path = "" + gpu_profile_requested = command_line%argument_present(["--gpu-profile"]) + if (args%order == 2 .and. ((allocated(dump_input_path) .and. len(dump_input_path) /= 0) .or. & + (allocated(args%dump_final_raw_path) .and. len(args%dump_final_raw_path) /= 0))) then + error stop "raw dumps require --order 4" + end if scalar_2D_initializer => scalar_field velocity_2D_initializer => taylor_green_velocity - s = scalar_2D_t(scalar_2D_initializer, order=4, cells=[51,51], x_min=[-pi,-pi], x_max=[pi,pi]) + call system_clock(count_rate=clock_rate) + if (clock_rate <= 0_int64) error stop "system_clock returned non-positive rate" - associate(v => vector_2D_t(velocity_2D_initializer, mold=s)) - - associate( & - scalar_file => s%to_file("scalar") & - ,velocity_file => v%to_file("velocity") & + call system_clock(count=setup_start) + if (allocated(args%input_path)) then + if (.not. args%no_csv) error stop "--input currently requires --no-csv" + call fa_read_raw_input(args%input_path, workspace_fields, args%cells(1), args%cells(2), args%steps, args%dt) + if (args%order /= default_order) error stop "--input requires order 4" + else + s = scalar_2D_t( & + scalar_2D_initializer & + ,order = args%order & + ,cells = args%cells & + ,x_min = x_min & + ,x_max = x_max & ) - call scalar_file%write_lines("example/scripts/scalar-initial.csv") - call velocity_file%write_lines("example/scripts/velocity.csv") - end associate + v = vector_2D_t(velocity_2D_initializer, mold=s) + s_initial = s + call fa_init_from_arrays(workspace_fields, & + real(s_initial%values(), kind=real32), & + real(v%values(x_dir), kind=real32), & + real(v%values(y_dir), kind=real32)) + end if + workspace_initial = workspace_fields%base + call system_clock(count=setup_end) - advance_time: & - block - real :: dt = 1E-4 - integer step - - do step = 1, 500 - associate(k1 => d_dt(s, v)) - associate(k2 => d_dt(s + (dt/2)*k1, v)) - associate(k3 => d_dt(s + (dt/2)*k2, v)) - associate(k4 => d_dt(s + dt*k3, v)) - s = s + (dt/6)*(k1 + 2*k2 + 2*k3 + k4) - end associate - end associate - end associate - end associate - end do - end block advance_time + initial_all_finite = all(ieee_is_finite(real(workspace_initial, kind=real64))) + if (.not. initial_all_finite) error stop "initial scalar field contains non-finite values" + if (.not. allocated(args%input_path)) then + call system_clock(count=initial_io_start) + if (allocated(dump_input_path) .and. len(dump_input_path) /= 0) call write_native_input(dump_input_path, s_initial, v, args) + if (.not. args%no_csv) then + associate(scalar_file => s_initial%to_file("scalar")) + call scalar_file%write_lines("example/scripts/scalar-initial.csv") + end associate + associate(velocity_file => v%to_file("velocity")) + call velocity_file%write_lines("example/scripts/velocity.csv") + end associate + end if + call system_clock(count=initial_io_end) + end if - associate(scalar_file => s%to_file("scalar")) - call scalar_file%write_lines("example/scripts/scalar-final.csv") - end associate + warmup_jit_seconds = 0E0 + if (gpu_profile_requested) then + call gpu_benchmark_begin(warmup_gpu_stats) + if (.not. warmup_gpu_stats%gpu_timing_available) error stop "gpu profiling requested but unavailable" + end if + if (args%solver == solver_workspace) then + call fa_advance_rk4(workspace_fields, args%cells(1), args%cells(2), 1, args%dt, pi, 0.5E0) + else + warmup_rhs = s + call advance_one_step(warmup_rhs, v, args%dt) + end if + if (args%solver == solver_workspace) then + warmup_checksum = sum(real(workspace_fields%base, kind=real64)) + else + warmup_checksum = sum(real(warmup_rhs%values(), kind=real64)) + end if + if (.not. ieee_is_finite(warmup_checksum)) error stop "warmup step contains non-finite values" + if (gpu_profile_requested) then + call gpu_benchmark_end(warmup_gpu_stats) + warmup_jit_seconds = warmup_gpu_stats%jit_seconds + end if + workspace_fields%base = workspace_initial + + if (args%solver == solver_workspace) then + if (gpu_profile_requested) then + call gpu_benchmark_begin(solve_gpu_stats) + if (.not. solve_gpu_stats%gpu_timing_available) error stop "gpu profiling requested but unavailable" + end if + call system_clock(count=solve_start) + call fa_advance_rk4(workspace_fields, args%cells(1), args%cells(2), args%steps, args%dt, pi, 0.5E0) + if (gpu_profile_requested) call gpu_benchmark_end(solve_gpu_stats) + call system_clock(count=solve_end) + else + s = s_initial + if (gpu_profile_requested) then + call gpu_benchmark_begin(solve_gpu_stats) + if (.not. solve_gpu_stats%gpu_timing_available) error stop "gpu profiling requested but unavailable" + end if + call system_clock(count=solve_start) + do step = 1, args%steps + call advance_one_step(s, v, args%dt) + end do + if (gpu_profile_requested) call gpu_benchmark_end(solve_gpu_stats) + call system_clock(count=solve_end) + workspace_fields%base = real(s%values(), kind=real32) + s_final = s + end if - end associate + final_all_finite = all(ieee_is_finite(real(workspace_fields%base, kind=real64))) + if (.not. final_all_finite) error stop "final scalar field contains non-finite values" + initial_checksum = sum(real(workspace_initial, kind=real64)) + final_checksum = sum(real(workspace_fields%base, kind=real64)) + l_inf_change = maxval(abs(real(workspace_fields%base, kind=real64) - real(workspace_initial, kind=real64))) + initial_norm_squared = sum(real(workspace_initial, kind=real64)**2) + difference_norm_squared = sum((real(workspace_fields%base, kind=real64) - real(workspace_initial, kind=real64))**2) + if (initial_norm_squared > 0E0) then + relative_l2_change = sqrt(difference_norm_squared / initial_norm_squared) + else + if (difference_norm_squared /= 0E0) error stop "zero scalar field changed" + relative_l2_change = 0E0 + end if + + if (.not. allocated(args%input_path)) then + call system_clock(count=final_io_start) + if (.not. args%no_csv) then + if (args%solver == solver_workspace) then + s_final = fa_reconstruct_scalar(workspace_fields%base, args%order, args%cells, x_min, x_max) + else + s_final = s + end if + associate(scalar_file => s_final%to_file("scalar")) + call scalar_file%write_lines("example/scripts/scalar-final.csv") + end associate + end if + call system_clock(count=final_io_end) + end if + + if (allocated(args%dump_final_raw_path) .and. len(args%dump_final_raw_path) /= 0) then + call fa_write_final_raw(args%dump_final_raw_path, workspace_fields, args%cells(1), args%cells(2), args%steps, args%dt) + end if + + setup_seconds = real(setup_end - setup_start, kind=real64) / real(clock_rate, kind=real64) + if (.not. allocated(args%input_path)) then + initial_io_seconds = real(initial_io_end - initial_io_start, kind=real64) / real(clock_rate, kind=real64) + final_io_seconds = real(final_io_end - final_io_start, kind=real64) / real(clock_rate, kind=real64) + else + initial_io_seconds = 0E0 + final_io_seconds = 0E0 + end if + solve_seconds = real(solve_end - solve_start, kind=real64) / real(clock_rate, kind=real64) + + write(output_unit,'(a)') "configuration:" + write(output_unit,'(a,i0)') " order = ", args%order + write(output_unit,'(a,i0,1x,i0)') " cells = ", args%cells(1), args%cells(2) + write(output_unit,'(a,i0)') " steps = ", args%steps + write(output_unit,'(a,es16.8)') " dt = ", args%dt + write(output_unit,'(a,1x,a)') "solver =", merge("workspace","operators", args%solver == solver_workspace) + write(output_unit,'(a,1x,g0)') "setup_seconds =", setup_seconds + write(output_unit,'(a,1x,g0)') "initial_io_seconds =", initial_io_seconds + write(output_unit,'(a,1x,g0)') "solve_seconds =", solve_seconds + write(output_unit,'(a,1x,g0)') "final_io_seconds =", final_io_seconds + write(output_unit,'(a,1x,l1)') "initial_all_finite =", all(ieee_is_finite(real(workspace_initial, kind=real64))) + write(output_unit,'(a,1x,l1)') "final_all_finite =", final_all_finite + write(output_unit,'(a,1x,g0)') "initial_checksum =", initial_checksum + write(output_unit,'(a,1x,g0)') "final_checksum =", final_checksum + write(output_unit,'(a,1x,g0)') "field_linf_change =", l_inf_change + write(output_unit,'(a,1x,g0)') "field_relative_l2_change =", relative_l2_change + if (allocated(dump_input_path) .and. len(dump_input_path) /= 0) write(output_unit,'(a,1x,a)') "dump_input_path =", dump_input_path + if (allocated(args%dump_final_raw_path) .and. len(args%dump_final_raw_path) /= 0) write(output_unit,'(a,1x,a)') "dump_final_raw_path =", args%dump_final_raw_path + if (gpu_profile_requested) then + write(output_unit,'(a,1x,g0)') "warmup_jit_seconds =", warmup_jit_seconds + write(output_unit,'(a,1x,l1)') "gpu_timing_available =", solve_gpu_stats%gpu_timing_available + write(output_unit,'(a,1x,i0)') "gpu_completed_command_buffers =", solve_gpu_stats%completed_command_buffers + write(output_unit,'(a,1x,i0)') "gpu_commits =", solve_gpu_stats%commits + write(output_unit,'(a,1x,i0)') "gpu_dispatches =", solve_gpu_stats%dispatches + write(output_unit,'(a,1x,i0)') "gpu_host_syncs =", solve_gpu_stats%host_syncs + write(output_unit,'(a,1x,i0)') "gpu_host_waits =", solve_gpu_stats%host_waits + write(output_unit,'(a,1x,i0)') "gpu_memcpy_h2d_bytes =", solve_gpu_stats%memcpy_h2d_bytes + write(output_unit,'(a,1x,i0)') "gpu_memcpy_d2h_bytes =", solve_gpu_stats%memcpy_d2h_bytes + write(output_unit,'(a,1x,i0)') "gpu_buffer_allocations =", solve_gpu_stats%buffer_allocations + write(output_unit,'(a,1x,i0)') "gpu_buffer_alloc_bytes =", solve_gpu_stats%buffer_alloc_bytes + write(output_unit,'(a,1x,g0)') "gpu_elapsed_seconds =", solve_gpu_stats%gpu_elapsed_seconds + write(output_unit,'(a,1x,g0)') "gpu_busy_seconds =", solve_gpu_stats%gpu_busy_seconds + write(output_unit,'(a,1x,g0)') "gpu_jit_seconds =", solve_gpu_stats%jit_seconds + end if contains @@ -92,4 +300,208 @@ pure function d_dt(s, v) result(ds_dt) ds_dt = .div. (D * .grad. s) - .div. (v * s) end function + subroutine advance_one_step(state, velocity, dt) + type(scalar_2D_t), intent(inout) :: state + type(vector_2D_t), intent(in) :: velocity + real, intent(in) :: dt + associate(k1 => d_dt(state, velocity)) + associate(k2 => d_dt(state + (dt/2)*k1, velocity)) + associate(k3 => d_dt(state + (dt/2)*k2, velocity)) + associate(k4 => d_dt(state + dt*k3, velocity)) + state = state + (dt/6)*(k1 + 2*k2 + 2*k3 + k4) + end associate + end associate + end associate + end associate + end subroutine + + subroutine reconstruct_final_scalar(values, final_scalar) + real(real32), intent(in) :: values(:,:) + type(scalar_2D_t), intent(out) :: final_scalar + final_scalar = fa_reconstruct_scalar(values, args%order, args%cells, x_min, x_max) + end subroutine + + function get_numerical_arguments() result(numerical_arguments) + type(numerical_arguments_t) numerical_arguments + character(len=:), allocatable :: cells_string, steps_string, dt_string, order_string + character(len=:), allocatable :: solver_string, input_string, dump_final_raw_string + + numerical_arguments%solver = 0 + numerical_arguments%benchmark = command_line%argument_present(["--benchmark"]) + numerical_arguments%no_csv = command_line%argument_present(["--no-csv"]) + if (numerical_arguments%benchmark) then + numerical_arguments%cells = [257,257] + numerical_arguments%steps = 30000 + end if + + cells_string = command_line%flag_value("--cells") + steps_string = command_line%flag_value("--steps") + dt_string = command_line%flag_value("--dt") + order_string = command_line%flag_value("--order") + solver_string = command_line%flag_value("--solver") + input_string = command_line%flag_value("--input") + dump_final_raw_string = command_line%flag_value("--dump-final-raw") + + if (len(cells_string) /= 0) then + call parse_int_pair(cells_string, numerical_arguments%cells) + numerical_arguments%cells_supplied = .true. + end if + if (len(steps_string) /= 0) then + call parse_positive_integer(steps_string, numerical_arguments%steps, "--steps") + numerical_arguments%steps_supplied = .true. + end if + if (len(dt_string) /= 0) then + call parse_positive_real(dt_string, numerical_arguments%dt, "--dt") + numerical_arguments%dt_supplied = .true. + end if + if (len(order_string) /= 0) call parse_order(order_string, numerical_arguments%order) + if (len(solver_string) /= 0) numerical_arguments%solver = parse_solver_mode(solver_string) + if (len(input_string) /= 0) numerical_arguments%input_path = input_string + if (len(dump_final_raw_string) /= 0) numerical_arguments%dump_final_raw_path = dump_final_raw_string + + if (allocated(numerical_arguments%input_path)) then + if (numerical_arguments%benchmark .or. numerical_arguments%cells_supplied .or. & + numerical_arguments%steps_supplied .or. numerical_arguments%dt_supplied) then + error stop "--input is exclusive with --benchmark, --cells, --steps, and --dt" + end if + if (numerical_arguments%order /= default_order) error stop "--input requires --order 4" + if (numerical_arguments%solver == 0) numerical_arguments%solver = solver_workspace + if (numerical_arguments%solver /= solver_workspace) error stop "--input requires --solver workspace" + end if + + if (numerical_arguments%solver == 0) then + if (numerical_arguments%order == 4) then + numerical_arguments%solver = solver_workspace + else + numerical_arguments%solver = solver_operators + end if + end if + + if (numerical_arguments%order == 2 .and. numerical_arguments%solver == solver_workspace) then + error stop "--solver workspace is unsupported for --order 2" + end if + + call validate_numerical_arguments(numerical_arguments) + + if (.not. allocated(numerical_arguments%input_path) .and. .not. numerical_arguments%dt_supplied) then + numerical_arguments%dt = default_dt * min(real(51, kind=real64) / real(numerical_arguments%cells(1), kind=real64), & + real(51, kind=real64) / real(numerical_arguments%cells(2), kind=real64))**2 + end if + end function + + subroutine validate_numerical_arguments(numerical_arguments) + type(numerical_arguments_t), intent(in) :: numerical_arguments + if (.not. any([numerical_arguments%order == 2, numerical_arguments%order == 4])) then + error stop "unsupported order; expected 2 or 4" + end if + if (any(numerical_arguments%cells <= 0)) error stop "cells must be positive" + if (any(numerical_arguments%cells < 2*numerical_arguments%order + 1)) error stop "cells must satisfy cells >= 2*order + 1" + if (numerical_arguments%steps <= 0) error stop "steps must be positive" + if (.not. allocated(numerical_arguments%input_path)) then + if (.not. numerical_arguments%dt_supplied .and. numerical_arguments%dt <= 0E0) error stop "dt must be positive" + if (numerical_arguments%dt <= 0E0) error stop "dt must be positive" + end if + end subroutine + + integer function parse_solver_mode(raw) result(mode) + character(len=*), intent(in) :: raw + character(len=:), allocatable :: trimmed + + trimmed = trim(adjustl(raw)) + select case (trimmed) + case ("workspace") + mode = solver_workspace + case ("operators") + mode = solver_operators + case default + error stop "--solver must be workspace or operators" + end select + end function + + subroutine parse_order(raw, value) + character(len=*), intent(in) :: raw + integer, intent(out) :: value + call read_integer_flag(raw, value, "--order") + if (.not. any([value == 2, value == 4])) error stop "order must be 2 or 4" + end subroutine + + subroutine parse_positive_integer(raw, value, flag_name) + character(len=*), intent(in) :: raw + integer, intent(out) :: value + character(len=*), intent(in) :: flag_name + call read_integer_flag(raw, value, flag_name) + if (value <= 0) error stop flag_name // " must be positive" + end subroutine + + subroutine parse_positive_real(raw, value, flag_name) + character(len=*), intent(in) :: raw, flag_name + real, intent(out) :: value + call read_real_flag(raw, value, flag_name) + if (value <= 0E0) error stop flag_name // " must be positive" + end subroutine + + subroutine read_integer_flag(raw, value, flag_name) + character(len=*), intent(in) :: raw, flag_name + integer, intent(out) :: value + integer :: ios + character(len=256) :: iomsg + character(len=:), allocatable :: buffer + + buffer = trim(adjustl(raw)) + + read(buffer, *, iostat=ios, iomsg=iomsg) value + if (ios /= 0) error stop flag_name // " parse error: " // trim(iomsg) + end subroutine + + subroutine read_real_flag(raw, value, flag_name) + character(len=*), intent(in) :: raw, flag_name + real, intent(out) :: value + integer :: ios + character(len=256) :: iomsg + character(len=:), allocatable :: buffer + + buffer = trim(adjustl(raw)) + + read(buffer, *, iostat=ios, iomsg=iomsg) value + if (ios /= 0) error stop flag_name // " parse error: " // trim(iomsg) + if (.not. ieee_is_finite(value)) error stop flag_name // " must be finite" + end subroutine + + subroutine parse_int_pair(raw, values) + character(len=*), intent(in) :: raw + integer, intent(out) :: values(2) + character(len=:), allocatable :: trimmed + integer :: comma, ios + character(len=256) :: iomsg + + trimmed = trim(adjustl(raw)) + comma = index(trimmed, ',') + if (comma == 0) then + call read_integer_flag(trimmed, values(1), "--cells") + values(2) = values(1) + else + read(trimmed, *, iostat=ios, iomsg=iomsg) values + if (ios /= 0) error stop "--cells parse error: " // trim(iomsg) + end if + end subroutine + + subroutine write_native_input(path, scalar_state, velocity_state, args) + character(len=*), intent(in) :: path + type(scalar_2D_t), intent(in) :: scalar_state + type(vector_2D_t), intent(in) :: velocity_state + type(numerical_arguments_t), intent(in) :: args + integer :: unit + character(len=256) :: header + + open(newunit=unit, file=path, access="stream", form="unformatted", status="replace", action="write") + write(unit) "FADV2D1" // new_line("a") + write(header,'(i0,1x,i0,1x,i0,1x,es24.16)') args%cells(1), args%cells(2), args%steps, real(args%dt, kind=real64) + write(unit) trim(header) // new_line("a") + write(unit) "mode raw" // new_line("a") + write(unit) real(scalar_state%values(), kind=real32) + write(unit) real(velocity_state%values(x_dir), kind=real32) + write(unit) real(velocity_state%values(y_dir), kind=real32) + close(unit) + end subroutine + end program diff --git a/example/formal_advection_workspace.F90 b/example/formal_advection_workspace.F90 new file mode 100644 index 0000000..69a0510 --- /dev/null +++ b/example/formal_advection_workspace.F90 @@ -0,0 +1,684 @@ +module formal_advection_workspace + use iso_fortran_env, only: real32, real64, int64 + use tensors_2D_m, only: scalar_2D_t + implicit none + private + public :: fa_config_t, fa_options_t, fa_fields_t, fa_parse_args + public :: fa_init_fields, fa_advance_rk4, fa_write_csv, fa_reconstruct_scalar + public :: fa_validate_cfg, fa_validate_working_set, fa_read_raw_input + public :: fa_now_seconds + + type :: fa_config_t + integer :: nx = 51 + integer :: ny = 51 + integer :: steps = 500 + real(real32) :: dt = 1.0e-4_real32 + end type + + type :: fa_options_t + type(fa_config_t) :: cfg + logical :: dump_csv = .false. + logical :: gpu_profile = .false. + logical :: dt_supplied = .false. + character(len=:), allocatable :: input_path + character(len=:), allocatable :: output_path + end type + + type :: fa_fields_t + real(real32), allocatable :: base(:,:) + real(real32), allocatable :: q(:,:) + real(real32), allocatable :: acc(:,:) + real(real32), allocatable :: vx(:,:) + real(real32), allocatable :: vy(:,:) + real(real32), allocatable :: diff_x(:,:) + real(real32), allocatable :: adv_x(:,:) + real(real32), allocatable :: diff_y(:,:) + real(real32), allocatable :: adv_y(:,:) + end type + + real(real32), allocatable, save :: fa_scalar_reconstruction_values(:,:) + +contains + + real(real64) function fa_now_seconds() result(now) + integer(int64) :: count, rate + call system_clock(count=count, count_rate=rate) + if (rate <= 0) error stop "system_clock is unavailable" + now = real(count, real64) / real(rate, real64) + end function fa_now_seconds + + pure real(real32) function fa_center_coord(i, n, minv, maxv) result(x) + integer, intent(in) :: i, n + real(real32), intent(in) :: minv, maxv + real(real32) :: dx + dx = (maxv - minv) / real(n, real32) + if (i == 1) then + x = minv + else if (i == n + 2) then + x = maxv + else + x = minv + (real(i, real32) - 1.5_real32) * dx + end if + end function fa_center_coord + + pure real(real32) function fa_face_coord(i, n, minv, maxv) result(x) + integer, intent(in) :: i, n + real(real32), intent(in) :: minv, maxv + real(real32) :: dx + dx = (maxv - minv) / real(n, real32) + x = minv + real(i - 1, real32) * dx + end function fa_face_coord + + pure real(real32) function fa_velocity_x(x, y) result(v) + real(real32), intent(in) :: x, y + v = 10.0_real32 * sin(x) * cos(y) + end function fa_velocity_x + + pure real(real32) function fa_velocity_y(x, y) result(v) + real(real32), intent(in) :: x, y + v = -10.0_real32 * cos(x) * sin(y) + end function fa_velocity_y + + subroutine fa_die(message) + character(len=*), intent(in) :: message + write(*, '(a)') trim(message) + error stop trim(message) + end subroutine fa_die + + integer(int64) function fa_checked_mul(a, b, what) result(product) + integer(int64), intent(in) :: a, b + character(len=*), intent(in) :: what + if (a /= 0_int64 .and. b > huge(0_int64) / a) then + call fa_die(trim(what) // " size overflow") + end if + product = a * b + end function fa_checked_mul + + subroutine fa_validate_working_set(nx, ny) + integer, intent(in) :: nx, ny + integer(int64) :: scalar_len, xface_len, yface_len, total_bytes + scalar_len = fa_checked_mul(int(nx + 2, int64), int(ny + 2, int64), & + "scalar field") + xface_len = fa_checked_mul(int(nx + 1, int64), int(ny + 2, int64), & + "x-face field") + yface_len = fa_checked_mul(int(nx + 2, int64), int(ny + 1, int64), & + "y-face field") + total_bytes = 3_int64 * scalar_len * 4_int64 & + + 3_int64 * xface_len * 4_int64 & + + 3_int64 * yface_len * 4_int64 + if (total_bytes > 512_int64 * 1024_int64 * 1024_int64) then + call fa_die("estimated working set exceeds 512 MiB") + end if + end subroutine fa_validate_working_set + + subroutine fa_validate_cfg(cfg, dt_supplied) + type(fa_config_t), intent(inout) :: cfg + logical, intent(in) :: dt_supplied + real(real32) :: scale + if (cfg%nx < 9 .or. cfg%ny < 9) then + call fa_die("cells must be at least 9 in each dimension") + end if + if (cfg%steps < 1) then + call fa_die("steps must be positive") + end if + if (.not. (cfg%dt > 0.0_real32)) then + call fa_die("dt must be finite and positive") + end if + if (.not. dt_supplied) then + scale = 51.0_real32 / real(max(cfg%nx, cfg%ny), real32) + cfg%dt = 1.0e-4_real32 * scale * scale + end if + end subroutine fa_validate_cfg + + pure real(real32) function fa_scalar_init(x, y, x0, y0, sigma) result(s) + real(real32), intent(in) :: x, y, x0, y0, sigma + real(real32) :: dx, dy, r + dx = x - x0 + dy = y - y0 + r = sqrt(dx * dx + dy * dy) + s = exp(-((r * r) / (2.0_real32 * sigma * sigma))) + end function fa_scalar_init + + subroutine fa_fill_analytic(fields, nx, ny, pi, x0, y0, sigma) + type(fa_fields_t), intent(inout) :: fields + integer, intent(in) :: nx, ny + real(real32), intent(in) :: pi, x0, y0, sigma + integer :: i, j + do concurrent (j = 1:ny + 2, i = 1:nx + 2) + fields%base(i, j) = fa_scalar_init( & + fa_center_coord(i, nx, -pi, pi), & + fa_center_coord(j, ny, -pi, pi), x0, y0, sigma) + fields%q(i, j) = fields%base(i, j) + fields%acc(i, j) = 0.0_real32 + end do + do concurrent (j = 1:ny + 2, i = 1:nx + 1) + fields%vx(i, j) = fa_velocity_x( & + fa_face_coord(i, nx, -pi, pi), & + fa_center_coord(j, ny, -pi, pi)) + end do + do concurrent (j = 1:ny + 1, i = 1:nx + 2) + fields%vy(i, j) = fa_velocity_y( & + fa_center_coord(i, nx, -pi, pi), & + fa_face_coord(j, ny, -pi, pi)) + end do + end subroutine fa_fill_analytic + + subroutine fa_read_line(unit, line, ios) + integer, intent(in) :: unit + character(len=:), allocatable, intent(out) :: line + integer, intent(out) :: ios + character(len=1) :: ch + character(len=1024) :: buffer + integer :: len_line + + buffer = '' + len_line = 0 + ios = 0 + do + read(unit, iostat=ios) ch + if (ios /= 0) exit + if (ch == new_line('a')) exit + len_line = len_line + 1 + if (len_line > len(buffer)) then + ios = 1 + exit + end if + buffer(len_line:len_line) = ch + end do + if (len_line > 0) then + line = buffer(:len_line) + else + line = '' + end if + end subroutine fa_read_line + + subroutine fa_read_raw_input(path, fields, nx, ny, steps, dt) + character(len=*), intent(in) :: path + type(fa_fields_t), intent(inout) :: fields + integer, intent(out) :: nx, ny, steps + real(real32), intent(out) :: dt + integer :: unit, ios + character(len=:), allocatable :: line + character(len=16) :: key, mode + character(len=512) :: msg + type(fa_config_t) :: cfg + + open(newunit=unit, file=path, status="old", action="read", & + access="stream", form="unformatted", iostat=ios, iomsg=msg) + if (ios /= 0) call fa_die(trim(msg)) + + call fa_read_line(unit, line, ios) + if (ios /= 0 .or. trim(line) /= "FADV2D1") then + call fa_die("invalid input header") + end if + call fa_read_line(unit, line, ios) + if (ios /= 0) call fa_die("missing input dimensions") + read(line, *, iostat=ios) cfg%nx, cfg%ny, cfg%steps, cfg%dt + if (ios /= 0) call fa_die("cannot parse input dimensions") + call fa_read_line(unit, line, ios) + if (ios /= 0) call fa_die("missing input mode") + read(line, *, iostat=ios) key, mode + if (ios /= 0 .or. trim(key) /= "mode") then + call fa_die("cannot parse input mode") + end if + if (trim(mode) /= "raw") then + call fa_die("unknown input mode") + end if + + call fa_validate_cfg(cfg, .true.) + nx = cfg%nx + ny = cfg%ny + steps = cfg%steps + dt = cfg%dt + call fa_validate_working_set(nx, ny) + call fa_allocate_fields(fields, nx, ny) + read(unit, iostat=ios) fields%base + if (ios /= 0) call fa_die("cannot read raw base field") + read(unit, iostat=ios) fields%vx + if (ios /= 0) call fa_die("cannot read raw vx field") + read(unit, iostat=ios) fields%vy + if (ios /= 0) call fa_die("cannot read raw vy field") + fields%q = fields%base + fields%acc = 0.0_real32 + close(unit) + end subroutine fa_read_raw_input + + pure real(real32) function fa_grad_face_x(q, nx, ny, i, j, inv_dx) result(g) + real(real32), intent(in) :: q(:,:) + integer, intent(in) :: nx, ny, i, j + real(real32), intent(in) :: inv_dx + integer :: n + n = nx + 2 + if (i == 1) then + g = inv_dx * ((-352.0_real32 / 105.0_real32) * q(1, j) & + + (35.0_real32 / 8.0_real32) * q(2, j) & + + (-35.0_real32 / 24.0_real32) * q(3, j) & + + (21.0_real32 / 40.0_real32) * q(4, j) & + + (-5.0_real32 / 56.0_real32) * q(5, j)) + else if (i == 2) then + g = inv_dx * ((16.0_real32 / 105.0_real32) * q(1, j) & + + (-31.0_real32 / 24.0_real32) * q(2, j) & + + (29.0_real32 / 24.0_real32) * q(3, j) & + + (-3.0_real32 / 40.0_real32) * q(4, j) & + + (1.0_real32 / 168.0_real32) * q(5, j)) + else if (i + 2 < n) then + g = inv_dx * ((1.0_real32 / 24.0_real32) * q(i - 1, j) & + + (-9.0_real32 / 8.0_real32) * q(i, j) & + + (9.0_real32 / 8.0_real32) * q(i + 1, j) & + + (-1.0_real32 / 24.0_real32) * q(i + 2, j)) + else if (i == n - 2) then + g = inv_dx * ((-1.0_real32 / 168.0_real32) * q(n - 4, j) & + + (3.0_real32 / 40.0_real32) * q(n - 3, j) & + + (-29.0_real32 / 24.0_real32) * q(n - 2, j) & + + (31.0_real32 / 24.0_real32) * q(n - 1, j) & + + (-16.0_real32 / 105.0_real32) * q(n, j)) + else + g = inv_dx * ((5.0_real32 / 56.0_real32) * q(n - 4, j) & + + (-21.0_real32 / 40.0_real32) * q(n - 3, j) & + + (35.0_real32 / 24.0_real32) * q(n - 2, j) & + + (-35.0_real32 / 8.0_real32) * q(n - 1, j) & + + (352.0_real32 / 105.0_real32) * q(n, j)) + end if + end function fa_grad_face_x + + pure real(real32) function fa_grad_face_y(q, nx, ny, i, j, inv_dy) result(g) + real(real32), intent(in) :: q(:,:) + integer, intent(in) :: nx, ny, i, j + real(real32), intent(in) :: inv_dy + integer :: n + n = ny + 2 + if (j == 1) then + g = inv_dy * ((-352.0_real32 / 105.0_real32) * q(i, 1) & + + (35.0_real32 / 8.0_real32) * q(i, 2) & + + (-35.0_real32 / 24.0_real32) * q(i, 3) & + + (21.0_real32 / 40.0_real32) * q(i, 4) & + + (-5.0_real32 / 56.0_real32) * q(i, 5)) + else if (j == 2) then + g = inv_dy * ((16.0_real32 / 105.0_real32) * q(i, 1) & + + (-31.0_real32 / 24.0_real32) * q(i, 2) & + + (29.0_real32 / 24.0_real32) * q(i, 3) & + + (-3.0_real32 / 40.0_real32) * q(i, 4) & + + (1.0_real32 / 168.0_real32) * q(i, 5)) + else if (j + 2 < n) then + g = inv_dy * ((1.0_real32 / 24.0_real32) * q(i, j - 1) & + + (-9.0_real32 / 8.0_real32) * q(i, j) & + + (9.0_real32 / 8.0_real32) * q(i, j + 1) & + + (-1.0_real32 / 24.0_real32) * q(i, j + 2)) + else if (j == n - 2) then + g = inv_dy * ((-1.0_real32 / 168.0_real32) * q(i, n - 4) & + + (3.0_real32 / 40.0_real32) * q(i, n - 3) & + + (-29.0_real32 / 24.0_real32) * q(i, n - 2) & + + (31.0_real32 / 24.0_real32) * q(i, n - 1) & + + (-16.0_real32 / 105.0_real32) * q(i, n)) + else + g = inv_dy * ((5.0_real32 / 56.0_real32) * q(i, n - 4) & + + (-21.0_real32 / 40.0_real32) * q(i, n - 3) & + + (35.0_real32 / 24.0_real32) * q(i, n - 2) & + + (-35.0_real32 / 8.0_real32) * q(i, n - 1) & + + (352.0_real32 / 105.0_real32) * q(i, n)) + end if + end function fa_grad_face_y + + pure real(real32) function fa_interp_face_x(q, nx, ny, i, j) result(v) + real(real32), intent(in) :: q(:,:) + integer, intent(in) :: nx, ny, i, j + integer :: n + n = nx + 2 + if (i == 1) then + v = q(1, j) + else if (i == 2) then + v = (-16.0_real32 / 112.0_real32) * q(1, j) & + + (70.0_real32 / 112.0_real32) * q(2, j) & + + (70.0_real32 / 112.0_real32) * q(3, j) & + + (-14.0_real32 / 112.0_real32) * q(4, j) & + + (2.0_real32 / 112.0_real32) * q(5, j) + else if (i + 2 < n) then + v = (-7.0_real32 / 112.0_real32) * q(i - 1, j) & + + (63.0_real32 / 112.0_real32) * q(i, j) & + + (63.0_real32 / 112.0_real32) * q(i + 1, j) & + + (-7.0_real32 / 112.0_real32) * q(i + 2, j) + else if (i == n - 2) then + v = (2.0_real32 / 112.0_real32) * q(n - 4, j) & + + (-14.0_real32 / 112.0_real32) * q(n - 3, j) & + + (70.0_real32 / 112.0_real32) * q(n - 2, j) & + + (70.0_real32 / 112.0_real32) * q(n - 1, j) & + + (-16.0_real32 / 112.0_real32) * q(n, j) + else + v = q(n, j) + end if + end function fa_interp_face_x + + pure real(real32) function fa_interp_face_y(q, nx, ny, i, j) result(v) + real(real32), intent(in) :: q(:,:) + integer, intent(in) :: nx, ny, i, j + integer :: n + n = ny + 2 + if (j == 1) then + v = q(i, 1) + else if (j == 2) then + v = (-16.0_real32 / 112.0_real32) * q(i, 1) & + + (70.0_real32 / 112.0_real32) * q(i, 2) & + + (70.0_real32 / 112.0_real32) * q(i, 3) & + + (-14.0_real32 / 112.0_real32) * q(i, 4) & + + (2.0_real32 / 112.0_real32) * q(i, 5) + else if (j + 2 < n) then + v = (-7.0_real32 / 112.0_real32) * q(i, j - 1) & + + (63.0_real32 / 112.0_real32) * q(i, j) & + + (63.0_real32 / 112.0_real32) * q(i, j + 1) & + + (-7.0_real32 / 112.0_real32) * q(i, j + 2) + else if (j == n - 2) then + v = (2.0_real32 / 112.0_real32) * q(i, n - 4) & + + (-14.0_real32 / 112.0_real32) * q(i, n - 3) & + + (70.0_real32 / 112.0_real32) * q(i, n - 2) & + + (70.0_real32 / 112.0_real32) * q(i, n - 1) & + + (-16.0_real32 / 112.0_real32) * q(i, n) + else + v = q(i, n) + end if + end function fa_interp_face_y + + pure real(real32) function fa_div_face_x(flux, nx, ny, i, j, inv_dx) result(d) + real(real32), intent(in) :: flux(:,:) + integer, intent(in) :: nx, ny, i, j + real(real32), intent(in) :: inv_dx + integer :: n + n = nx + 1 + if (i == 1) then + d = inv_dx * ((-11.0_real32 / 12.0_real32) * flux(1, j) & + + (17.0_real32 / 24.0_real32) * flux(2, j) & + + (3.0_real32 / 8.0_real32) * flux(3, j) & + + (-5.0_real32 / 24.0_real32) * flux(4, j) & + + (1.0_real32 / 24.0_real32) * flux(5, j)) + else if (i + 1 < n) then + d = inv_dx * ((1.0_real32 / 24.0_real32) * flux(i - 1, j) & + + (-9.0_real32 / 8.0_real32) * flux(i, j) & + + (9.0_real32 / 8.0_real32) * flux(i + 1, j) & + + (-1.0_real32 / 24.0_real32) * flux(i + 2, j)) + else + d = inv_dx * ((-1.0_real32 / 24.0_real32) * flux(n - 4, j) & + + (5.0_real32 / 24.0_real32) * flux(n - 3, j) & + + (-3.0_real32 / 8.0_real32) * flux(n - 2, j) & + + (-17.0_real32 / 24.0_real32) * flux(n - 1, j) & + + (11.0_real32 / 12.0_real32) * flux(n, j)) + end if + end function fa_div_face_x + + pure real(real32) function fa_div_face_y(flux, nx, ny, i, j, inv_dy) result(d) + real(real32), intent(in) :: flux(:,:) + integer, intent(in) :: nx, ny, i, j + real(real32), intent(in) :: inv_dy + integer :: n + n = ny + 1 + if (j == 1) then + d = inv_dy * ((-11.0_real32 / 12.0_real32) * flux(i, 1) & + + (17.0_real32 / 24.0_real32) * flux(i, 2) & + + (3.0_real32 / 8.0_real32) * flux(i, 3) & + + (-5.0_real32 / 24.0_real32) * flux(i, 4) & + + (1.0_real32 / 24.0_real32) * flux(i, 5)) + else if (j + 1 < n) then + d = inv_dy * ((1.0_real32 / 24.0_real32) * flux(i, j - 1) & + + (-9.0_real32 / 8.0_real32) * flux(i, j) & + + (9.0_real32 / 8.0_real32) * flux(i, j + 1) & + + (-1.0_real32 / 24.0_real32) * flux(i, j + 2)) + else + d = inv_dy * ((-1.0_real32 / 24.0_real32) * flux(i, n - 4) & + + (5.0_real32 / 24.0_real32) * flux(i, n - 3) & + + (-3.0_real32 / 8.0_real32) * flux(i, n - 2) & + + (-17.0_real32 / 24.0_real32) * flux(i, n - 1) & + + (11.0_real32 / 12.0_real32) * flux(i, n)) + end if + end function fa_div_face_y + + subroutine fa_allocate_fields(fields, nx, ny) + type(fa_fields_t), intent(inout) :: fields + integer, intent(in) :: nx, ny + allocate(fields%base(1:nx+2, 1:ny+2)) + allocate(fields%q(1:nx+2, 1:ny+2)) + allocate(fields%acc(1:nx+2, 1:ny+2)) + allocate(fields%vx(1:nx+1, 1:ny+2)) + allocate(fields%vy(1:nx+2, 1:ny+1)) + allocate(fields%diff_x(1:nx+1, 1:ny+2)) + allocate(fields%adv_x(1:nx+1, 1:ny+2)) + allocate(fields%diff_y(1:nx+2, 1:ny+1)) + allocate(fields%adv_y(1:nx+2, 1:ny+1)) + end subroutine fa_allocate_fields + + subroutine fa_init_fields(fields, nx, ny, pi, x0, y0, sigma) + type(fa_fields_t), intent(inout) :: fields + integer, intent(in) :: nx, ny + real(real32), intent(in) :: pi, x0, y0, sigma + call fa_allocate_fields(fields, nx, ny) + call fa_fill_analytic(fields, nx, ny, pi, x0, y0, sigma) + end subroutine fa_init_fields + + subroutine fa_init_from_arrays(fields, base, vx, vy) + type(fa_fields_t), intent(inout) :: fields + real(real32), intent(in) :: base(:,:), vx(:,:), vy(:,:) + integer :: nx, ny + + nx = size(base, 1) - 2 + ny = size(base, 2) - 2 + call fa_allocate_fields(fields, nx, ny) + fields%base = base + fields%q = base + fields%acc = 0.0_real32 + fields%vx = vx + fields%vy = vy + fields%diff_x = 0.0_real32 + fields%adv_x = 0.0_real32 + fields%diff_y = 0.0_real32 + fields%adv_y = 0.0_real32 + end subroutine fa_init_from_arrays + + subroutine fa_write_final_raw(path, fields, nx, ny, steps, dt) + character(len=*), intent(in) :: path + type(fa_fields_t), intent(in) :: fields + integer, intent(in) :: nx, ny, steps + real(real32), intent(in) :: dt + integer :: unit + character(len=256) :: header + + open(newunit=unit, file=path, access="stream", form="unformatted", & + status="replace", action="write") + write(unit) "FADV2D1" // new_line("a") + write(header,'(i0,1x,i0,1x,i0,1x,es24.16)') nx, ny, steps, real(dt, real64) + write(unit) trim(header) // new_line("a") + write(unit) "mode raw-final" // new_line("a") + write(unit) real(fields%base, kind=real32) + close(unit) + end subroutine fa_write_final_raw + + pure function fa_scalar_reconstruction_initializer(x, y) result(field) + real, intent(in) :: x(:), y(:) + real field(size(x), size(y)) + + if (.not. allocated(fa_scalar_reconstruction_values)) then + error stop "scalar reconstruction values are unavailable" + end if + field = real(fa_scalar_reconstruction_values, kind=kind(field)) + end function fa_scalar_reconstruction_initializer + + function fa_reconstruct_scalar(values, order, cells, x_min, x_max) result(scalar) + real(real32), intent(in) :: values(:,:) + integer, intent(in) :: order + integer, intent(in) :: cells(:) + real(real32), intent(in) :: x_min(:), x_max(:) + type(scalar_2D_t) :: scalar + + fa_scalar_reconstruction_values = values + scalar = scalar_2D_t(fa_scalar_reconstruction_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) + deallocate(fa_scalar_reconstruction_values) + end function fa_reconstruct_scalar + + subroutine fa_advance_rk4(fields, nx, ny, steps, dt, pi, diffusion) + type(fa_fields_t), intent(inout) :: fields + integer, intent(in) :: nx, ny, steps + real(real32), intent(in) :: dt, pi, diffusion + integer :: step, stage, i, j + real(real32) :: inv_dx, inv_dy, half_dt, dt_over_6, rhs, total + real(real32), allocatable :: vx_local(:,:), vy_local(:,:) + + inv_dx = real(nx, real32) / (2.0_real32 * pi) + inv_dy = real(ny, real32) / (2.0_real32 * pi) + half_dt = 0.5_real32 * dt + dt_over_6 = dt / 6.0_real32 + allocate(vx_local(1:nx + 1, 1:ny + 2)) + allocate(vy_local(1:nx + 2, 1:ny + 1)) + do concurrent (j = 1:ny + 2, i = 1:nx + 1) + vx_local(i, j) = fields%vx(i, j) + end do + do concurrent (j = 1:ny + 1, i = 1:nx + 2) + vy_local(i, j) = fields%vy(i, j) + end do + do concurrent (j = 1:ny + 2, i = 1:nx + 2) + fields%q(i, j) = fields%base(i, j) + fields%acc(i, j) = 0.0_real32 + end do + + do step = 1, steps + do stage = 1, 4 + do concurrent (j = 1:ny + 2, i = 1:nx + 2) + if (i <= nx + 1) then + fields%diff_x(i, j) = diffusion * fa_grad_face_x( & + fields%q, nx, ny, i, j, inv_dx) + fields%adv_x(i, j) = vx_local(i, j) * fa_interp_face_x( & + fields%q, nx, ny, i, j) + end if + if (j <= ny + 1) then + fields%diff_y(i, j) = diffusion * fa_grad_face_y( & + fields%q, nx, ny, i, j, inv_dy) + fields%adv_y(i, j) = vy_local(i, j) * fa_interp_face_y( & + fields%q, nx, ny, i, j) + end if + end do + + do concurrent (j = 2:ny + 1, i = 2:nx + 1) local(rhs, total) + rhs = (fa_div_face_x(fields%diff_x, nx, ny, i - 1, j, & + inv_dx) + fa_div_face_y(fields%diff_y, nx, ny, i, j - 1, & + inv_dy)) - (fa_div_face_x(fields%adv_x, nx, ny, i - 1, j, & + inv_dx) + fa_div_face_y(fields%adv_y, nx, ny, i, j - 1, & + inv_dy)) + select case (stage) + case (1) + fields%acc(i, j) = rhs + fields%q(i, j) = fields%base(i, j) + half_dt * rhs + case (2) + fields%acc(i, j) = fields%acc(i, j) + 2.0_real32 * rhs + fields%q(i, j) = fields%base(i, j) + half_dt * rhs + case (3) + fields%acc(i, j) = fields%acc(i, j) + 2.0_real32 * rhs + fields%q(i, j) = fields%base(i, j) + dt * rhs + case default + total = fields%acc(i, j) + rhs + fields%base(i, j) = fields%base(i, j) + dt_over_6 * total + fields%q(i, j) = fields%base(i, j) + end select + end do + end do + end do + deallocate(vx_local) + deallocate(vy_local) + end subroutine fa_advance_rk4 + + subroutine fa_write_csv(path, fields, nx, ny) + character(len=*), intent(in) :: path + type(fa_fields_t), intent(in) :: fields + integer, intent(in) :: nx, ny + integer :: unit, i, j, ios + character(len=512) :: message + real(real32) :: pi + + pi = acos(-1.0_real32) + + open(newunit=unit, file=path, status="replace", action="write", & + form="formatted", iostat=ios, iomsg=message) + if (ios /= 0) error stop message + write(unit, '(a)') 'x,y,scalar' + do j = 1, ny + 2 + do i = 1, nx + 2 + write(unit, '(f0.9, ",", f0.9, ",", es0.8)') & + fa_center_coord(i, nx, -pi, pi), & + fa_center_coord(j, ny, -pi, pi), & + fields%base(i, j) + end do + write(unit, *) + end do + close(unit) + end subroutine fa_write_csv + + subroutine fa_parse_args(opts) + type(fa_options_t), intent(inout) :: opts + integer :: argc, i + character(len=256) :: arg, value + + argc = command_argument_count() + i = 1 + do while (i <= argc) + call get_command_argument(i, arg) + select case (trim(arg)) + case ("--cells") + i = i + 1 + call get_command_argument(i, value) + read(value, *) opts%cfg%nx, opts%cfg%ny + case ("--steps") + i = i + 1 + call get_command_argument(i, value) + read(value, *) opts%cfg%steps + case ("--dt") + i = i + 1 + call get_command_argument(i, value) + read(value, *) opts%cfg%dt + opts%dt_supplied = .true. + case ("--gpu-profile") + opts%gpu_profile = .true. + case ("--input") + i = i + 1 + call get_command_argument(i, value) + opts%input_path = trim(value) + case ("--dump", "--csv") + opts%dump_csv = .true. + if (i < argc) then + call get_command_argument(i + 1, value) + if (value(1:1) /= '-') then + opts%output_path = trim(value) + i = i + 1 + end if + end if + case ("--output") + i = i + 1 + call get_command_argument(i, value) + opts%output_path = trim(value) + opts%dump_csv = .true. + case ("--help", "-h") + write(*, '(a)') "usage: fortran_workspace [--cells nx,ny] " // & + "[--steps n] [--dt dt] [--input path] " // & + "[--dump [csv]] [--output path] [--gpu-profile]" + stop 0 + case default + if (arg(1:1) /= '-') then + ! Positional numeric convenience: nx ny steps dt. + if (i == 1 .and. argc >= 4) then + read(arg, *) opts%cfg%nx + call get_command_argument(2, value) + read(value, *) opts%cfg%ny + call get_command_argument(3, value) + read(value, *) opts%cfg%steps + call get_command_argument(4, value) + read(value, *) opts%cfg%dt + opts%dt_supplied = .true. + i = 4 + else + error stop "unknown argument" + end if + else + error stop "unknown argument" + end if + end select + i = i + 1 + end do + if (.not. allocated(opts%output_path)) then + opts%output_path = "build/fortran/final.csv" + end if + end subroutine fa_parse_args + +end module formal_advection_workspace diff --git a/example/gpu_benchmark_m.F90 b/example/gpu_benchmark_m.F90 new file mode 100644 index 0000000..274057f --- /dev/null +++ b/example/gpu_benchmark_m.F90 @@ -0,0 +1,123 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +module gpu_benchmark_m + use iso_c_binding, only : c_double, c_int, c_int64_t + use iso_fortran_env, only : int64, real64 + implicit none + private + + type, public :: gpu_benchmark_stats_t + logical :: gpu_timing_available = .false. + integer(int64) :: completed_command_buffers = 0 + integer(int64) :: commits = 0 + integer(int64) :: dispatches = 0 + integer(int64) :: host_syncs = 0 + integer(int64) :: host_waits = 0 + integer(int64) :: memcpy_h2d_bytes = 0 + integer(int64) :: memcpy_d2h_bytes = 0 + integer(int64) :: buffer_allocations = 0 + integer(int64) :: buffer_alloc_bytes = 0 + real(real64) :: gpu_elapsed_seconds = 0E0 + real(real64) :: gpu_busy_seconds = 0E0 + real(real64) :: jit_seconds = 0E0 + end type + + public :: gpu_benchmark_begin, gpu_benchmark_end + + interface + subroutine formal_gpu_benchmark_begin(gpu_timing_available) bind(C, name="formal_gpu_benchmark_begin") + import c_int + integer(c_int), intent(out) :: gpu_timing_available + end subroutine + + subroutine formal_gpu_benchmark_end( & + gpu_timing_available & + ,completed_command_buffers & + ,commits & + ,dispatches & + ,host_syncs & + ,host_waits & + ,memcpy_h2d_bytes & + ,memcpy_d2h_bytes & + ,buffer_allocations & + ,buffer_alloc_bytes & + ,gpu_elapsed_seconds & + ,gpu_busy_seconds & + ,jit_seconds & + ) bind(C, name="formal_gpu_benchmark_end") + import c_double, c_int, c_int64_t + integer(c_int), intent(out) :: gpu_timing_available + integer(c_int64_t), intent(out) :: completed_command_buffers + integer(c_int64_t), intent(out) :: commits + integer(c_int64_t), intent(out) :: dispatches + integer(c_int64_t), intent(out) :: host_syncs + integer(c_int64_t), intent(out) :: host_waits + integer(c_int64_t), intent(out) :: memcpy_h2d_bytes + integer(c_int64_t), intent(out) :: memcpy_d2h_bytes + integer(c_int64_t), intent(out) :: buffer_allocations + integer(c_int64_t), intent(out) :: buffer_alloc_bytes + real(c_double), intent(out) :: gpu_elapsed_seconds + real(c_double), intent(out) :: gpu_busy_seconds + real(c_double), intent(out) :: jit_seconds + end subroutine + end interface + +contains + + subroutine gpu_benchmark_begin(stats) + type(gpu_benchmark_stats_t), intent(out) :: stats + integer(c_int) :: gpu_timing_available + + call formal_gpu_benchmark_begin(gpu_timing_available) + stats = gpu_benchmark_stats_t(gpu_timing_available = gpu_timing_available /= 0) + end subroutine + + subroutine gpu_benchmark_end(stats) + type(gpu_benchmark_stats_t), intent(out) :: stats + integer(c_int) :: gpu_timing_available + integer(c_int64_t) :: completed_command_buffers + integer(c_int64_t) :: commits + integer(c_int64_t) :: dispatches + integer(c_int64_t) :: host_syncs + integer(c_int64_t) :: host_waits + integer(c_int64_t) :: memcpy_h2d_bytes + integer(c_int64_t) :: memcpy_d2h_bytes + integer(c_int64_t) :: buffer_allocations + integer(c_int64_t) :: buffer_alloc_bytes + real(c_double) :: gpu_elapsed_seconds + real(c_double) :: gpu_busy_seconds + real(c_double) :: jit_seconds + + call formal_gpu_benchmark_end( & + gpu_timing_available & + ,completed_command_buffers & + ,commits & + ,dispatches & + ,host_syncs & + ,host_waits & + ,memcpy_h2d_bytes & + ,memcpy_d2h_bytes & + ,buffer_allocations & + ,buffer_alloc_bytes & + ,gpu_elapsed_seconds & + ,gpu_busy_seconds & + ,jit_seconds & + ) + + stats%gpu_timing_available = gpu_timing_available /= 0 + stats%completed_command_buffers = completed_command_buffers + stats%commits = commits + stats%dispatches = dispatches + stats%host_syncs = host_syncs + stats%host_waits = host_waits + stats%memcpy_h2d_bytes = memcpy_h2d_bytes + stats%memcpy_d2h_bytes = memcpy_d2h_bytes + stats%buffer_allocations = buffer_allocations + stats%buffer_alloc_bytes = buffer_alloc_bytes + stats%gpu_elapsed_seconds = gpu_elapsed_seconds + stats%gpu_busy_seconds = gpu_busy_seconds + stats%jit_seconds = jit_seconds + end subroutine + +end module gpu_benchmark_m diff --git a/example/gpu_benchmark_shim.c b/example/gpu_benchmark_shim.c new file mode 100644 index 0000000..ac9a457 --- /dev/null +++ b/example/gpu_benchmark_shim.c @@ -0,0 +1,158 @@ +#include +#include + +#ifdef FORMAL_ENABLE_GPU_BENCHMARK +#include "lfortran_gpu_runtime.h" + +static lfortran_gpu_ctx *gpu_ctx = NULL; +static lfortran_gpu_metal_telemetry begin_snapshot; +static int has_begin_snapshot = 0; + +static void ensure_context(void) +{ + if (gpu_ctx == NULL) { + gpu_ctx = lfortran_gpu_init(); + } +} + +void formal_gpu_benchmark_begin(int *gpu_timing_available) +{ + ensure_context(); + memset(&begin_snapshot, 0, sizeof(begin_snapshot)); + if (gpu_ctx != NULL) { + lfortran_gpu_metal_telemetry_snapshot_sync(gpu_ctx, &begin_snapshot); + } + has_begin_snapshot = 1; + if (gpu_timing_available != NULL) { + *gpu_timing_available = begin_snapshot.gpu_timing_available; + } +} + +void formal_gpu_benchmark_end( + int *gpu_timing_available, + long long *completed_command_buffers, + long long *commits, + long long *dispatches, + long long *host_syncs, + long long *host_waits, + long long *memcpy_h2d_bytes, + long long *memcpy_d2h_bytes, + long long *buffer_allocations, + long long *buffer_alloc_bytes, + double *gpu_elapsed_seconds, + double *gpu_busy_seconds, + double *jit_seconds) +{ + lfortran_gpu_metal_telemetry end_snapshot; + + ensure_context(); + memset(&end_snapshot, 0, sizeof(end_snapshot)); + if (gpu_ctx != NULL) { + lfortran_gpu_metal_telemetry_snapshot_sync(gpu_ctx, &end_snapshot); + } + + if (gpu_timing_available != NULL) { + *gpu_timing_available = end_snapshot.gpu_timing_available; + } + if (completed_command_buffers != NULL) { + *completed_command_buffers = (long long)(end_snapshot.completed_command_buffers - begin_snapshot.completed_command_buffers); + } + if (commits != NULL) { + *commits = (long long)(end_snapshot.commits - begin_snapshot.commits); + } + if (dispatches != NULL) { + *dispatches = (long long)(end_snapshot.dispatches - begin_snapshot.dispatches); + } + if (host_syncs != NULL) { + *host_syncs = (long long)(end_snapshot.host_syncs - begin_snapshot.host_syncs); + } + if (host_waits != NULL) { + *host_waits = (long long)(end_snapshot.host_waits - begin_snapshot.host_waits); + } + if (memcpy_h2d_bytes != NULL) { + *memcpy_h2d_bytes = (long long)(end_snapshot.memcpy_h2d_bytes - begin_snapshot.memcpy_h2d_bytes); + } + if (memcpy_d2h_bytes != NULL) { + *memcpy_d2h_bytes = (long long)(end_snapshot.memcpy_d2h_bytes - begin_snapshot.memcpy_d2h_bytes); + } + if (buffer_allocations != NULL) { + *buffer_allocations = (long long)(end_snapshot.buffer_allocations - begin_snapshot.buffer_allocations); + } + if (buffer_alloc_bytes != NULL) { + *buffer_alloc_bytes = (long long)(end_snapshot.buffer_alloc_bytes - begin_snapshot.buffer_alloc_bytes); + } + if (gpu_elapsed_seconds != NULL) { + *gpu_elapsed_seconds = end_snapshot.gpu_elapsed_seconds - begin_snapshot.gpu_elapsed_seconds; + } + if (gpu_busy_seconds != NULL) { + *gpu_busy_seconds = end_snapshot.gpu_busy_seconds - begin_snapshot.gpu_busy_seconds; + } + if (jit_seconds != NULL) { + *jit_seconds = end_snapshot.jit_seconds - begin_snapshot.jit_seconds; + } + has_begin_snapshot = 0; +} +#else +void formal_gpu_benchmark_begin(int *gpu_timing_available) +{ + if (gpu_timing_available != NULL) { + *gpu_timing_available = 0; + } +} + +void formal_gpu_benchmark_end( + int *gpu_timing_available, + long long *completed_command_buffers, + long long *commits, + long long *dispatches, + long long *host_syncs, + long long *host_waits, + long long *memcpy_h2d_bytes, + long long *memcpy_d2h_bytes, + long long *buffer_allocations, + long long *buffer_alloc_bytes, + double *gpu_elapsed_seconds, + double *gpu_busy_seconds, + double *jit_seconds) +{ + if (gpu_timing_available != NULL) { + *gpu_timing_available = 0; + } + if (completed_command_buffers != NULL) { + *completed_command_buffers = 0; + } + if (commits != NULL) { + *commits = 0; + } + if (dispatches != NULL) { + *dispatches = 0; + } + if (host_syncs != NULL) { + *host_syncs = 0; + } + if (host_waits != NULL) { + *host_waits = 0; + } + if (memcpy_h2d_bytes != NULL) { + *memcpy_h2d_bytes = 0; + } + if (memcpy_d2h_bytes != NULL) { + *memcpy_d2h_bytes = 0; + } + if (buffer_allocations != NULL) { + *buffer_allocations = 0; + } + if (buffer_alloc_bytes != NULL) { + *buffer_alloc_bytes = 0; + } + if (gpu_elapsed_seconds != NULL) { + *gpu_elapsed_seconds = 0.0; + } + if (gpu_busy_seconds != NULL) { + *gpu_busy_seconds = 0.0; + } + if (jit_seconds != NULL) { + *jit_seconds = 0.0; + } +} +#endif diff --git a/example/scripts/compare_2d_fields.py b/example/scripts/compare_2d_fields.py new file mode 100755 index 0000000..39e6253 --- /dev/null +++ b/example/scripts/compare_2d_fields.py @@ -0,0 +1,221 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import csv +import math +import shutil +import sys +from pathlib import Path + + +def finite_float(raw: str, path: Path, line_number: int, kind: str) -> float: + try: + value = float(raw) + except ValueError as exc: + raise SystemExit(f"{path}:{line_number}: invalid {kind}: {raw!r}") from exc + if not math.isfinite(value): + raise SystemExit(f"{path}:{line_number}: non-finite {kind}") + return value + + +def load_field(path: Path): + if not path.exists(): + raise SystemExit(f"{path}: missing file") + + xs = [] + ys = [] + x_seen = set() + y_seen = set() + values = {} + with path.open(newline="") as handle: + reader = csv.reader(handle) + header = next(reader, None) + if header is None: + raise SystemExit(f"{path}: empty file") + if len(header) < 3: + raise SystemExit(f"{path}: expected x, y, value columns") + for line_number, row in enumerate(reader, start=2): + if not row: + continue + if len(row) != 3: + raise SystemExit(f"{path}:{line_number}: expected exactly 3 CSV columns") + x = finite_float(row[0], path, line_number, "x") + y = finite_float(row[1], path, line_number, "y") + value = finite_float(row[2], path, line_number, "value") + if x not in x_seen: + xs.append(x) + x_seen.add(x) + if y not in y_seen: + ys.append(y) + y_seen.add(y) + key = (x, y) + if key in values: + raise SystemExit(f"{path}:{line_number}: duplicate grid point {x},{y}") + values[key] = value + + if not xs or not ys: + raise SystemExit(f"{path}: empty field data") + xs.sort() + ys.sort() + if len(values) != len(xs) * len(ys): + raise SystemExit(f"{path}: incomplete field data") + + field = [[values[(x, y)] for y in ys] for x in xs] + return xs, ys, field + + +def compare_fields(reference, observed, rtol: float, atol: float): + if len(reference) != len(observed) or len(reference[0]) != len(observed[0]): + raise SystemExit("field shapes differ") + + max_abs_diff = 0.0 + max_rel_diff = 0.0 + changed = False + for i in range(len(reference)): + for j in range(len(reference[0])): + ref = reference[i][j] + obs = observed[i][j] + diff = abs(obs - ref) + scale = max(abs(ref), abs(obs)) + max_abs_diff = max(max_abs_diff, diff) + max_rel_diff = max(max_rel_diff, diff / scale if scale else 0.0) + changed = changed or diff != 0.0 + if diff > (atol + rtol * scale): + raise SystemExit( + f"field mismatch at ({i},{j}): abs diff {diff:.16e} exceeds " + f"atol+rtol*scale = {(atol + rtol * scale):.16e}" + ) + ref_sq = sum(reference[i][j] * reference[i][j] for i in range(len(reference)) for j in range(len(reference[0]))) + diff_sq = sum((observed[i][j] - reference[i][j]) ** 2 for i in range(len(reference)) for j in range(len(reference[0]))) + relative_l2 = math.sqrt(diff_sq / ref_sq) if ref_sq else (0.0 if diff_sq == 0.0 else math.inf) + return max_abs_diff, max_rel_diff, relative_l2, changed + + +def summarize_field_differences(reference, observed): + if len(reference) != len(observed) or len(reference[0]) != len(observed[0]): + raise SystemExit("field shapes differ") + + max_abs_diff = 0.0 + max_rel_diff = 0.0 + changed = False + for i in range(len(reference)): + for j in range(len(reference[0])): + ref = reference[i][j] + obs = observed[i][j] + diff = abs(obs - ref) + scale = max(abs(ref), abs(obs)) + max_abs_diff = max(max_abs_diff, diff) + max_rel_diff = max(max_rel_diff, diff / scale if scale else 0.0) + changed = changed or diff != 0.0 + ref_sq = sum(reference[i][j] * reference[i][j] for i in range(len(reference)) for j in range(len(reference[0]))) + diff_sq = sum((observed[i][j] - reference[i][j]) ** 2 for i in range(len(reference)) for j in range(len(reference[0]))) + relative_l2 = math.sqrt(diff_sq / ref_sq) if ref_sq else (0.0 if diff_sq == 0.0 else math.inf) + return max_abs_diff, max_rel_diff, relative_l2, changed + + +def load_and_validate(reference_path: Path, observed_path: Path): + reference_x, reference_y, reference = load_field(reference_path) + observed_x, observed_y, observed = load_field(observed_path) + if reference_x != observed_x or reference_y != observed_y: + raise SystemExit("coordinate grids differ") + return reference, observed + + +def run_comparison(reference_path: Path, observed_path: Path, rtol: float, atol: float, evolution: bool): + reference, observed = load_and_validate(reference_path, observed_path) + if evolution: + max_abs_diff, max_rel_diff, relative_l2, changed = summarize_field_differences(reference, observed) + else: + max_abs_diff, max_rel_diff, relative_l2, changed = compare_fields(reference, observed, rtol, atol) + print(f"max_abs_diff={max_abs_diff:.16e}") + print(f"max_rel_diff={max_rel_diff:.16e}") + print(f"relative_l2={relative_l2:.16e}") + print(f"changed={str(changed).lower()}") + if evolution and not changed: + raise SystemExit("field did not change") + + +def write_csv(path: Path, xs, ys, func): + path.parent.mkdir(parents=True, exist_ok=True) + with path.open("w", newline="") as handle: + writer = csv.writer(handle) + writer.writerow(["x", "y", "value"]) + for x in xs: + for y in ys: + writer.writerow([f"{x:.16e}", f"{y:.16e}", f"{func(x, y):.16e}"]) + + +def expect_failure(label: str, func): + try: + func() + except SystemExit: + return + raise SystemExit(f"self-test failed: expected failure for {label}") + + +def self_test(): + repo_root = Path(__file__).resolve().parents[2] + scratch = repo_root / "build" / "formal-gpu" / "compare_2d_fields-selftest" + if scratch.exists(): + shutil.rmtree(scratch) + scratch.mkdir(parents=True) + + xs = [0.0, 1.0] + ys = [0.0, 2.0] + reference = scratch / "reference.csv" + candidate = scratch / "candidate.csv" + shifted = scratch / "shifted.csv" + shifted_grid = scratch / "shifted-grid.csv" + malformed = scratch / "malformed.csv" + nonfinite = scratch / "nonfinite.csv" + duplicate = scratch / "duplicate.csv" + empty = scratch / "empty.csv" + + write_csv(reference, xs, ys, lambda x, y: 0.0) + write_csv(candidate, xs, ys, lambda x, y: 0.0) + write_csv(shifted, xs, ys, lambda x, y: 1.0 if x == 1.0 and y == 2.0 else 0.0) + write_csv(shifted_grid, [0.0, 1.1], ys, lambda x, y: 0.0) + + malformed.write_text("x,y,value\n0,0\n", encoding="utf-8") + nonfinite.write_text("x,y,value\n0,0,nan\n", encoding="utf-8") + duplicate.write_text("x,y,value\n0,0,1\n0,0,2\n", encoding="utf-8") + empty.write_text("", encoding="utf-8") + + run_comparison(reference, candidate, 0.0, 0.0, False) + run_comparison(reference, shifted, 0.0, 0.0, True) + expect_failure("mismatched values", lambda: run_comparison(reference, shifted, 0.0, 0.0, False)) + expect_failure("mismatched grid", lambda: run_comparison(reference, shifted_grid, 0.0, 0.0, False)) + expect_failure("malformed csv", lambda: load_field(malformed)) + expect_failure("nonfinite csv", lambda: load_field(nonfinite)) + expect_failure("duplicate csv", lambda: load_field(duplicate)) + expect_failure("empty csv", lambda: load_field(empty)) + print("self-test ok") + + +def main(argv): + parser = argparse.ArgumentParser(description="Compare 2D field CSVs") + parser.add_argument("reference", nargs="?") + parser.add_argument("candidate", nargs="?") + parser.add_argument("--rtol", type=float, default=0.0) + parser.add_argument("--atol", type=float, default=0.0) + parser.add_argument("--evolution", action="store_true", help="require candidate to differ from reference") + parser.add_argument("--self-test", action="store_true", help="run comparator self-tests") + args = parser.parse_args(argv[1:]) + + if args.self_test: + self_test() + return 0 + + if args.reference is None or args.candidate is None: + parser.error("reference and candidate CSV paths are required unless --self-test is used") + + if not math.isfinite(args.rtol) or not math.isfinite(args.atol) or args.rtol < 0.0 or args.atol < 0.0: + raise SystemExit("--rtol and --atol must be finite and non-negative") + + run_comparison(Path(args.reference), Path(args.candidate), args.rtol, args.atol, args.evolution) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main(sys.argv)) diff --git a/example/scripts/validate_2d_advection_workspace.py b/example/scripts/validate_2d_advection_workspace.py new file mode 100644 index 0000000..ea824a5 --- /dev/null +++ b/example/scripts/validate_2d_advection_workspace.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import struct +import subprocess +import sys +from pathlib import Path + +from compare_2d_fields import run_comparison as compare_csv + + +def parse_col_major_payload(payload: bytes, nx: int, ny: int): + count = (nx + 2) * (ny + 2) + if len(payload) != count * 4: + raise SystemExit(f"unexpected raw payload size: expected {count * 4} bytes") + values = struct.unpack(f"<{count}f", payload) + return [[values[i + j * (nx + 2)] for j in range(ny + 2)] for i in range(nx + 2)] + + +def parse_raw_headered(path: Path): + with path.open("rb") as handle: + header = handle.readline().decode("ascii").strip() + if header != "FADV2D1": + raise SystemExit(f"{path}: invalid header") + dims = handle.readline().decode("ascii").strip().split() + if len(dims) != 4: + raise SystemExit(f"{path}: invalid dimensions line") + nx, ny, steps = map(int, dims[:3]) + dt = float(dims[3]) + mode = handle.readline().decode("ascii").strip() + if mode not in {"mode raw-final", "mode raw"}: + raise SystemExit(f"{path}: invalid mode {mode!r}") + payload = handle.read() + return nx, ny, steps, dt, parse_col_major_payload(payload, nx, ny) + + +def parse_raw_headerless(path: Path, nx: int, ny: int): + payload = path.read_bytes() + return parse_col_major_payload(payload, nx, ny) + + +def compare_fields(reference, observed, rtol: float, atol: float): + if len(reference) != len(observed) or len(reference[0]) != len(observed[0]): + raise SystemExit("raw field shapes differ") + max_abs = 0.0 + max_rel = 0.0 + changed = False + for i in range(len(reference)): + for j in range(len(reference[0])): + ref = reference[i][j] + obs = observed[i][j] + diff = abs(obs - ref) + scale = max(abs(ref), abs(obs)) + max_abs = max(max_abs, diff) + max_rel = max(max_rel, diff / scale if scale else 0.0) + changed = changed or diff != 0.0 + if diff > atol + rtol * scale: + raise SystemExit( + f"raw field mismatch at ({i},{j}): abs diff {diff:.16e} exceeds " + f"atol+rtol*scale = {(atol + rtol * scale):.16e}" + ) + print(f"raw_max_abs_diff={max_abs:.16e}") + print(f"raw_max_rel_diff={max_rel:.16e}") + print(f"raw_changed={str(changed).lower()}") + + +def run(command, cwd: Path): + subprocess.run(command, cwd=cwd, check=True) + + +def main(argv): + parser = argparse.ArgumentParser(description="Validate the 2D advection-diffusion example") + parser.add_argument("--executable", required=True) + parser.add_argument("--gpu-runner", nargs="+", required=True, help="command prefix that owns the executable") + parser.add_argument("--workdir", default="build/formal-gpu/workspace-validation") + parser.add_argument("--cells", default="9,11") + parser.add_argument("--steps", default="3") + parser.add_argument("--dt", default="1e-4") + parser.add_argument("--rtol", type=float, default=1e-6) + parser.add_argument("--atol", type=float, default=1e-8) + parser.add_argument("--raw-input") + parser.add_argument("--reference-raw") + parser.add_argument("--native-final-raw") + parser.add_argument("--gpu-profile", action="store_true") + args = parser.parse_args(argv[1:]) + + if args.rtol < 0.0 or args.atol < 0.0: + raise SystemExit("--rtol and --atol must be non-negative") + + workdir = Path(args.workdir) + workdir.mkdir(parents=True, exist_ok=True) + exe = Path(args.executable).resolve() + repo_root = Path(__file__).resolve().parents[2] + runner = list(args.gpu_runner) + + if args.raw_input: + final_raw = workdir / "workspace.final.raw" + command = runner + [str(exe), "--solver", "workspace", "--input", args.raw_input, "--no-csv", "--dump-final-raw", str(final_raw)] + if args.gpu_profile: + command.append("--gpu-profile") + run(command, repo_root) + candidate_nx, candidate_ny, candidate_steps, candidate_dt, candidate = parse_raw_headered(final_raw) + if args.native_final_raw: + native = parse_raw_headerless(Path(args.native_final_raw), candidate_nx, candidate_ny) + compare_fields(candidate, native, args.rtol, args.atol) + if args.reference_raw: + _, _, _, _, reference = parse_raw_headered(Path(args.reference_raw)) + compare_fields(reference, candidate, args.rtol, args.atol) + return 0 + + operators_prefix = workdir / "operators" + workspace_prefix = workdir / "workspace" + operators_final = operators_prefix.with_suffix(".scalar-final.csv") + workspace_final = workspace_prefix.with_suffix(".scalar-final.csv") + operators_initial = operators_prefix.with_suffix(".scalar-initial.csv") + workspace_initial = workspace_prefix.with_suffix(".scalar-initial.csv") + + for prefix in (operators_prefix, workspace_prefix): + for suffix in ("scalar-initial.csv", "scalar-final.csv", "velocity.csv"): + p = prefix.parent / f"{prefix.name}.{suffix}" + if p.exists(): + p.unlink() + + operators_cmd = runner + [str(exe), "--solver", "operators", "--cells", args.cells, "--steps", args.steps, "--dt", args.dt] + workspace_cmd = runner + [str(exe), "--solver", "workspace", "--cells", args.cells, "--steps", args.steps, "--dt", args.dt] + if args.gpu_profile: + operators_cmd.append("--gpu-profile") + workspace_cmd.append("--gpu-profile") + run(operators_cmd, repo_root) + (workdir / "operators.scalar-initial.csv").write_bytes((repo_root / "example/scripts/scalar-initial.csv").read_bytes()) + (workdir / "operators.scalar-final.csv").write_bytes((repo_root / "example/scripts/scalar-final.csv").read_bytes()) + run(workspace_cmd, repo_root) + (workdir / "workspace.scalar-initial.csv").write_bytes((repo_root / "example/scripts/scalar-initial.csv").read_bytes()) + (workdir / "workspace.scalar-final.csv").write_bytes((repo_root / "example/scripts/scalar-final.csv").read_bytes()) + compare_csv(operators_final, workspace_final, rtol=args.rtol, atol=args.atol, evolution=False) + compare_csv(operators_initial, workspace_initial, rtol=args.rtol, atol=args.atol, evolution=False) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main(sys.argv)) From 197d18c96d91a8de31a6ca76fb93dc4fcdb678e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 20:21:46 -0600 Subject: [PATCH 05/18] Add native Metal advection benchmark (cherry picked from commit 48e70be3fbb4591789c55f249b4fa7da44a23955) --- README.md | 23 +- benchmarks/2d-advection-metal/.gitignore | 6 + .../2d-advection-metal/LFortran-NOTICE.txt | 45 + benchmarks/2d-advection-metal/PERFORMANCE.md | 67 + benchmarks/2d-advection-metal/README.md | 85 ++ .../2d-advection-metal/formal_advection.m | 1200 +++++++++++++++++ .../2d-advection-metal/formal_advection.metal | 254 ++++ benchmarks/2d-advection-metal/run.sh | 69 + benchmarks/2d-advection-metal/run_formal.sh | 50 + 9 files changed, 1790 insertions(+), 9 deletions(-) create mode 100644 benchmarks/2d-advection-metal/.gitignore create mode 100644 benchmarks/2d-advection-metal/LFortran-NOTICE.txt create mode 100644 benchmarks/2d-advection-metal/PERFORMANCE.md create mode 100644 benchmarks/2d-advection-metal/README.md create mode 100644 benchmarks/2d-advection-metal/formal_advection.m create mode 100644 benchmarks/2d-advection-metal/formal_advection.metal create mode 100755 benchmarks/2d-advection-metal/run.sh create mode 100755 benchmarks/2d-advection-metal/run_formal.sh diff --git a/README.md b/README.md index e748482..d74a94d 100644 --- a/README.md +++ b/README.md @@ -82,17 +82,22 @@ For GPU benchmarking, compile with `-DFORMAL_ENABLE_GPU_BENCHMARK`, pass `--gpu-profile`, and set `LFORTRAN_GPU_PROFILE=1` to collect Metal telemetry through the runtime snapshot shim. -Build LFortran with Pixi, then run its application benchmark entry point. -It supplies the C-shim flags, checks the actual Metal executable, and runs -it under the GPU guard: +The native C + Metal reference benchmark is kept outside FPM's automatic +example discovery under +[`benchmarks/2d-advection-metal`](./benchmarks/2d-advection-metal). It does +not add an Objective-C `main` to the ordinary Formal examples or tests. +After building LFortran with Pixi, run the native reference or the actual +Formal GPU example through its wrappers: ```bash -export FORMAL_ROOT=/path/to/formal -cd /path/to/lfortran -pixi run -e llvm11 bash examples/formal_advection/run_formal.sh +export LFORTRAN_ROOT=/path/to/lfortran +bash benchmarks/2d-advection-metal/run.sh +# or: +bash benchmarks/2d-advection-metal/run_formal.sh ``` -The default is the measured 2048 x 2048, 1000-step case; append example -arguments to select another workload. Build directories are keyed by the -compiler binary hash to avoid stale kernels after compiler changes. +The native wrapper defaults to a CPU smoke case; Metal is explicitly guarded +and the measured 2048 x 2048, 1000-step case is selected by passing its +arguments. See the benchmark README for the exact provenance, source guard, +and reproducibility contract. ### Other example programs See this repository's [example](./example) subdirectory for additional demonstrations of diff --git a/benchmarks/2d-advection-metal/.gitignore b/benchmarks/2d-advection-metal/.gitignore new file mode 100644 index 0000000..9208bbf --- /dev/null +++ b/benchmarks/2d-advection-metal/.gitignore @@ -0,0 +1,6 @@ +/build/ +*.bin +*.csv +*.dSYM/ +*.fadv2d +*.log diff --git a/benchmarks/2d-advection-metal/LFortran-NOTICE.txt b/benchmarks/2d-advection-metal/LFortran-NOTICE.txt new file mode 100644 index 0000000..352dfa7 --- /dev/null +++ b/benchmarks/2d-advection-metal/LFortran-NOTICE.txt @@ -0,0 +1,45 @@ +The files formal_advection.m and formal_advection.metal are derived from +the LFortran repository, commit 3d4c7941c5a153adb996ed62309225e871c06d9b, +under the following BSD 3-Clause terms. + +Copyright (c) 2019-2020, Triad National Security, LLC. All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this +list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the following disclaimer in the documentation +and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its contributors +may be used to endorse or promote products derived from this software without +specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. + +The LFortran source repository also contains the following government-rights +notice: + +This program was produced under U.S. Government contract +89233218CNA000001 for Los Alamos National Laboratory (LANL), which is +operated by Triad National Security, LLC for the U.S. Department of +Energy/National Nuclear Security Administration. All rights in the program +are reserved by Triad National Security, LLC, and the U.S. Department of +Energy/National Nuclear Security Administration. The Government is granted +for itself and others acting on its behalf a nonexclusive, paid-up, +irrevocable worldwide license in this material to reproduce, prepare +derivative works, distribute copies to the public, perform publicly and +display publicly, and to permit others to do so. diff --git a/benchmarks/2d-advection-metal/PERFORMANCE.md b/benchmarks/2d-advection-metal/PERFORMANCE.md new file mode 100644 index 0000000..5aad44f --- /dev/null +++ b/benchmarks/2d-advection-metal/PERFORMANCE.md @@ -0,0 +1,67 @@ +# Native C + Metal performance notes + +The native program is a reproducible oracle and a performance reference for +Formal's corrected 2D advection-diffusion workspace. It uses FP32 mimetic +gradient, interpolation, and divergence operators with classical RK4. It +does not replace the method with a five-point Laplacian. + +## Verified scheduling + +The Metal host keeps the verified bounded schedule: + +- four RK4 stages per timestep; +- eight timesteps and 64 dispatches per command buffer; +- at most two command buffers in flight; +- an autorelease pool per encoded command buffer; +- serial retirement and unioned command-buffer timestamps for `gpu_busy`. + +The CPU implementation is intentionally retained as the oracle. The +corrected one-sided gradient rows, fixed extended scalar boundaries, and +`k1 + 2*k2 + 2*k3 + k4` accumulation order are shared by the CPU and Metal +implementations. + +## Reproduction + +On the Apple M4 Max used for the original comparison, the measured Gaussian +case was 2048 x 2048 cells for 1000 steps with the example's scaled timestep +(approximately `6.20126741e-8`). Run the native reference with: + +```bash +LFORTRAN_ROOT=/path/to/lfortran \ + bash benchmarks/2d-advection-metal/run.sh \ + --mode metal --cells 2048,2048 --steps 1000 +``` + +The wrapper compiles with the configured macOS SDK compiler and invokes +LFortran's existing GPU guard directly. The guard is required for every +Metal run and serializes users through `/tmp/gpurun.lock.d`. +The native program writes no CSV unless `--dump` or `--output` is supplied; +use `--dump-bin` when a full-precision field is needed for comparison. + +The original paired measurement used the same raw Formal input for both +implementations: + +| variant | solve wall | GPU command-buffer time | +|---|---:|---:| +| Formal + LFortran | 2.195701 s | 2.186008 s | +| native C + Metal | 2.250356 s | 2.247697 s | + +Those values are historical verification samples, not a promise of a fixed +runtime across machines. The native result was within the measured Formal +throughput, with largest elementwise difference +`5.96046448e-8` and relative L2 difference `5.98075916e-9`. + +## Index experiment + +`--index-mode linear` preserves the signed flattened-index implementation +from the verified source. On the smaller 257 x 257 cache-resident workload +it measured about 1.55 s GPU versus 1.08 s for direct 2D indexing. An +unsigned-index experiment was slightly slower and was rejected; this document +does not claim that unsigned indexing is faster. + +## Provenance + +The native source and shader are from LFortran commit +`3d4c7941c5a153adb996ed62309225e871c06d9b`; see the benchmark +[`README.md`](./README.md) and [`LFortran-NOTICE.txt`](./LFortran-NOTICE.txt) +for the exact source paths and license handling. diff --git a/benchmarks/2d-advection-metal/README.md b/benchmarks/2d-advection-metal/README.md new file mode 100644 index 0000000..4a9832e --- /dev/null +++ b/benchmarks/2d-advection-metal/README.md @@ -0,0 +1,85 @@ +# Native 2D advection-diffusion reference + +This directory contains the native C + Metal reference implementation used +to measure the Formal 2D advection-diffusion example. It is deliberately +under `benchmarks/`, outside FPM's automatic `example/` discovery, so the +Objective-C `main` is never added to ordinary Formal examples or tests. + +The CPU path is the numerical oracle. Both paths implement the corrected +order-4 mimetic gradient, center-to-face interpolation, and divergence +closures followed by classical RK4: + +```text +rhs = div(D * grad(s)) - div(v * interp(s)) +``` + +This is not a five-point Laplacian or a lower-order substitute. Scalar +boundaries remain fixed, the FP32 field layout matches Formal, and the RK4 +accumulator is updated in the same `k1 + 2*k2 + 2*k3 + k4` order. + +## Build and run + +The wrappers use the LFortran `llvm11` Pixi environment. They accept an +explicit `LFORTRAN_ROOT`; otherwise they look for a sibling checkout at +`../lfortran` and fail with an actionable message if it is absent: + +```bash +LFORTRAN_ROOT=/path/to/lfortran \ + bash benchmarks/2d-advection-metal/run.sh +``` + +The wrapper re-enters Pixi with +`pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11`, selects the +configured macOS SDK compiler with `xcrun`, and builds +`build/formal_advection`. With no arguments it runs a small CPU smoke case. +The native executable reads `formal_advection.metal` from its current working +directory; `run.sh` changes to this directory before invoking it, so it cannot +accidentally load a stale shader. + +Metal is opt-in and doubly guarded: + +```bash +FA_ALLOW_METAL=1 LFORTRAN_ROOT=/path/to/lfortran \ + bash benchmarks/2d-advection-metal/run.sh \ + --mode metal --cells 257,257 --steps 30000 +``` + +Metal runs invoke the existing LFortran +`examples/nbody/gpu_run.sh` directly. That preserves its shared +`/tmp/gpurun.lock.d` serialization, timeout, and memory guard; there is no +unguarded fallback. `GPU_TIMEOUT` and `GPU_MAXMEM_MB` are passed through to +that guard. + +Useful native options include `--dump `, `--dump-bin [path]`, +`--index-mode 2d|linear`, and `--tg x[,y]`. The `linear` mode is the signed +flattened-index experiment. The verified comparison found it slower on the +small cache-resident case; no claim that unsigned indexing is faster is made. + +To run the actual Formal GPU example rather than the native oracle: + +```bash +LFORTRAN_ROOT=/path/to/lfortran \ + bash benchmarks/2d-advection-metal/run_formal.sh +``` + +That wrapper compiles the existing +`example/2D-advection-diffusion.F90` with LFortran, checks Metal linkage, and +uses the same guarded launch path. It does not add native sources to FPM. + +## Provenance and licensing + +`formal_advection.m` and `formal_advection.metal` were copied verbatim from +the LFortran repository at commit +`3d4c7941c5a153adb996ed62309225e871c06d9b` +(`examples/formal_advection/formal_advection.m` and +`formal_advection.metal`). That commit is the source of the verified CPU +oracle, Metal kernels, two-deep scheduling, and signed-index experiment. +The adapted wrappers and this documentation are the Formal publication +integration. The upstream LFortran BSD terms are retained in +[`LFortran-NOTICE.txt`](./LFortran-NOTICE.txt); the surrounding Formal tree +remains covered by [`LICENSE.txt`](../../LICENSE.txt). + +Generated binaries, logs, CSV/raw fields, and build caches belong under +`build/` and are ignored. No Fortran workspace fixture is duplicated here: +the Formal implementation remains the tracked +`example/2D-advection-diffusion.F90`. diff --git a/benchmarks/2d-advection-metal/formal_advection.m b/benchmarks/2d-advection-metal/formal_advection.m new file mode 100644 index 0000000..a2d79a6 --- /dev/null +++ b/benchmarks/2d-advection-metal/formal_advection.m @@ -0,0 +1,1200 @@ +#import +#import + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +static const float FA_PI = 3.14159265358979323846f; +static const float FA_DIFFUSION = 0.5f; +static const float FA_SIGMA = 3.14159265358979323846f / 8.0f; +static const float FA_X0 = -3.14159265358979323846f / 2.0f; +static const float FA_Y0 = -3.14159265358979323846f / 2.0f; +static const char *FA_FORMAL_51_COORDS[] = { + "-3.141592741013", "-3.079992771149", "-2.956793069839", "-2.833593368530", + "-2.710393667221", "-2.587193965912", "-2.463994264603", "-2.340794563293", + "-2.217594623566", "-2.094395160675", "-1.971195340157", "-1.847995638847", + "-1.724795937538", "-1.601596117020", "-1.478396415710", "-1.355196714401", + "-1.231997013092", "-1.108797311783", "-0.9855976104736", "-0.8623979091644", + "-0.7391982078552", "-0.6159985065460", "-0.4927988052368", "-0.3695991039276", + "-0.2463994026184", "-0.1231994628906", "0.2384185791016E-06", "0.1231999397278", + "0.2463996410370", "0.3695993423462", "0.4927990436554", "0.6159987449646", + "0.7391984462738", "0.8623981475830", "0.9855980873108", "1.108797550201", + "1.231997489929", "1.355196952820", "1.478396892548", "1.601596355438", + "1.724796295166", "1.847995758057", "1.971195697784", "2.094395160675", + "2.217595100403", "2.340794563293", "2.463994503021", "2.587193965912", + "2.710393905640", "2.833593845367", "2.956793308258", "3.079993247986", + "3.141592741013", +}; + +typedef struct { + uint32_t nx; + uint32_t ny; + uint32_t steps; + float dt; +} fa_config_t; + +typedef struct { + size_t scalar_len; + size_t xface_len; + size_t yface_len; + float *base; + float *stage[2]; + float *accum; + float *vx; + float *vy; +} fa_fields_t; + +typedef struct { + bool metal; + bool dump_csv; + bool dump_bin; + char *input_path; + char *output_path; + char *binary_output_path; + fa_config_t cfg; + bool dt_supplied; + uint32_t tg_x; + uint32_t tg_y; + bool tg_supplied; + bool index_mode_supplied; + enum { + FA_INDEX_DIRECT_2D, + FA_INDEX_LINEAR, + } index_mode; +} fa_options_t; + +typedef struct { + double jit_seconds; + double solve_seconds; + double gpu_busy; + bool final_all_finite; +} fa_run_stats_t; + +static double fa_now_seconds(void) +{ + static mach_timebase_info_data_t timebase = {0, 0}; + if (timebase.denom == 0) { + if (mach_timebase_info(&timebase) != KERN_SUCCESS) { + fprintf(stderr, "cannot query mach timebase\n"); + exit(1); + } + } + return (double)mach_absolute_time() * (double)timebase.numer / + (double)timebase.denom * 1e-9; +} + +static void fa_die(const char *message) +{ + fprintf(stderr, "%s\n", message); + exit(1); +} + +static void fa_dief(const char *fmt, ...) +{ + va_list ap; + va_start(ap, fmt); + vfprintf(stderr, fmt, ap); + fprintf(stderr, "\n"); + va_end(ap); + exit(1); +} + +static size_t fa_checked_mul_size(size_t a, size_t b, const char *what) +{ + if (a != 0 && b > SIZE_MAX / a) { + fa_dief("%s size overflow", what); + } + return a * b; +} + +static const char *fa_center_coord_text(size_t i, uint32_t n) +{ + if (n == 51 && i < sizeof(FA_FORMAL_51_COORDS) / sizeof(FA_FORMAL_51_COORDS[0])) { + return FA_FORMAL_51_COORDS[i]; + } + return NULL; +} + +static float fa_center_coord(size_t i, uint32_t n, float min, float max) +{ + const char *text = fa_center_coord_text(i, n); + if (text) { + return strtof(text, NULL); + } + volatile float dx = (max - min) / (float)n; + if (i == 0) { + return min; + } + if (i == (size_t)n + 1) { + return max; + } + volatile float base = min + 0.5f * dx; + volatile float offset = ((float)i - 1.0f) * dx; + return base + offset; +} + +static float fa_face_coord(size_t i, uint32_t n, float min, float max) +{ + volatile float dx = (max - min) / (float)n; + if (i == 0) { + return min; + } + if (i == (size_t)n) { + return max; + } + return min + (float)i * dx; +} + +static size_t fa_scalar_index(size_t i, size_t j, uint32_t nx) +{ + return i + (size_t)(nx + 2) * j; +} + +static size_t fa_xface_index(size_t i, size_t j, uint32_t nx) +{ + return i + (size_t)(nx + 1) * j; +} + +static size_t fa_yface_index(size_t i, size_t j, uint32_t nx) +{ + return i + (size_t)(nx + 2) * j; +} + +static float fa_scalar_init(float x, float y) +{ + const float dx = x - FA_X0; + const float dy = y - FA_Y0; + const float r = sqrtf(dx * dx + dy * dy); + return expf(-((r * r) / (2.0f * FA_SIGMA * FA_SIGMA))); +} + +static float fa_velocity_x(float x, float y) +{ + return 10.0f * sinf(x) * cosf(y); +} + +static float fa_velocity_y(float x, float y) +{ + return -10.0f * cosf(x) * sinf(y); +} + +static void fa_fields_free(fa_fields_t *fields) +{ + free(fields->base); + free(fields->stage[0]); + free(fields->stage[1]); + free(fields->accum); + free(fields->vx); + free(fields->vy); + memset(fields, 0, sizeof(*fields)); +} + +static void fa_fields_alloc(fa_fields_t *fields, const fa_config_t *cfg) +{ + memset(fields, 0, sizeof(*fields)); + fields->scalar_len = fa_checked_mul_size((size_t)cfg->nx + 2, + (size_t)cfg->ny + 2, "scalar field"); + fields->xface_len = fa_checked_mul_size((size_t)cfg->nx + 1, + (size_t)cfg->ny + 2, "x-face field"); + fields->yface_len = fa_checked_mul_size((size_t)cfg->nx + 2, + (size_t)cfg->ny + 1, "y-face field"); + + fields->base = calloc(fields->scalar_len, sizeof(float)); + fields->stage[0] = calloc(fields->scalar_len, sizeof(float)); + fields->stage[1] = calloc(fields->scalar_len, sizeof(float)); + fields->accum = calloc(fields->scalar_len, sizeof(float)); + fields->vx = calloc(fields->xface_len, sizeof(float)); + fields->vy = calloc(fields->yface_len, sizeof(float)); + if (!fields->base || !fields->stage[0] || !fields->stage[1] || + !fields->accum || !fields->vx || !fields->vy) { + fa_die("cannot allocate fields"); + } +} + +static void fa_generate_analytic(const fa_config_t *cfg, fa_fields_t *fields) +{ + for (size_t j = 0; j < (size_t)cfg->ny + 2; ++j) { + const float y = fa_center_coord(j, cfg->ny, -FA_PI, FA_PI); + for (size_t i = 0; i < (size_t)cfg->nx + 2; ++i) { + const float x = fa_center_coord(i, cfg->nx, -FA_PI, FA_PI); + const size_t idx = fa_scalar_index(i, j, cfg->nx); + fields->base[idx] = fa_scalar_init(x, y); + fields->stage[0][idx] = fields->base[idx]; + fields->stage[1][idx] = fields->base[idx]; + } + } + + for (size_t j = 0; j < (size_t)cfg->ny + 2; ++j) { + const float y = fa_center_coord(j, cfg->ny, -FA_PI, FA_PI); + for (size_t i = 0; i < (size_t)cfg->nx + 1; ++i) { + const float x = fa_face_coord(i, cfg->nx, -FA_PI, FA_PI); + fields->vx[fa_xface_index(i, j, cfg->nx)] = fa_velocity_x(x, y); + } + } + + for (size_t j = 0; j < (size_t)cfg->ny + 1; ++j) { + const float y = fa_face_coord(j, cfg->ny, -FA_PI, FA_PI); + for (size_t i = 0; i < (size_t)cfg->nx + 2; ++i) { + const float x = fa_center_coord(i, cfg->nx, -FA_PI, FA_PI); + fields->vy[fa_yface_index(i, j, cfg->nx)] = fa_velocity_y(x, y); + } + } +} + +static void fa_validate_cfg(fa_config_t *cfg, bool dt_supplied) +{ + if (cfg->nx < 9 || cfg->ny < 9) { + fa_die("cells must be at least 9 in each dimension"); + } + if (cfg->steps == 0) { + fa_die("steps must be positive"); + } + if (!(isfinite(cfg->dt) && cfg->dt > 0.0f)) { + fa_die("dt must be finite and positive"); + } + if (!dt_supplied) { + const float sx = 51.0f / (float)cfg->nx; + const float sy = 51.0f / (float)cfg->ny; + const float scale = fminf(sx, sy); + cfg->dt = 1.0e-4f * scale * scale; + } +} + +static uint32_t fa_parse_u32_strict(const char *raw, const char *flag_name) +{ + errno = 0; + char *end = NULL; + unsigned long value = strtoul(raw, &end, 10); + if (errno != 0 || end == raw || *end != '\0' || value == 0 || + value > UINT32_MAX) { + fa_dief("%s parse error: expected positive unsigned integer", flag_name); + } + return (uint32_t)value; +} + +static void fa_validate_working_set(const fa_config_t *cfg) +{ + const size_t scalar_len = fa_checked_mul_size((size_t)cfg->nx + 2, + (size_t)cfg->ny + 2, "scalar field"); + const size_t xface_len = fa_checked_mul_size((size_t)cfg->nx + 1, + (size_t)cfg->ny + 2, "x-face field"); + const size_t yface_len = fa_checked_mul_size((size_t)cfg->nx + 2, + (size_t)cfg->ny + 1, "y-face field"); + const size_t total_bytes = 4 * scalar_len * sizeof(float) + + 2 * xface_len * sizeof(float) + + 2 * yface_len * sizeof(float); + if (total_bytes > 512u * 1024u * 1024u) { + fa_die("estimated GPU working set exceeds 512 MiB"); + } +} + +static float fa_parse_float_strict(const char *raw, const char *flag_name) +{ + errno = 0; + char *end = NULL; + float value = strtof(raw, &end); + if (errno != 0 || end == raw || *end != '\0' || !isfinite(value)) { + fa_dief("%s parse error: expected finite real", flag_name); + } + return value; +} + +static void fa_read_raw_input(const char *path, fa_config_t *cfg, fa_fields_t *fields) +{ + FILE *f = fopen(path, "rb"); + if (!f) { + fa_dief("cannot open input file: %s", path); + } + + char line[256]; + if (!fgets(line, sizeof(line), f) || strncmp(line, "FADV2D1", 7) != 0) { + fclose(f); + fa_die("invalid input header"); + } + if (!fgets(line, sizeof(line), f)) { + fclose(f); + fa_die("missing input dimensions"); + } + if (sscanf(line, "%u %u %u %f", &cfg->nx, &cfg->ny, &cfg->steps, &cfg->dt) != 4) { + fclose(f); + fa_die("cannot parse input dimensions"); + } + if (!fgets(line, sizeof(line), f)) { + fclose(f); + fa_die("missing input mode"); + } + char mode[64] = {0}; + if (sscanf(line, "mode %63s", mode) != 1) { + fclose(f); + fa_die("cannot parse input mode"); + } + fa_validate_cfg(cfg, true); + fa_validate_working_set(cfg); + fa_fields_alloc(fields, cfg); + if (strcmp(mode, "raw") == 0) { + if (fread(fields->base, sizeof(float), fields->scalar_len, f) != + fields->scalar_len || + fread(fields->vx, sizeof(float), fields->xface_len, f) != + fields->xface_len || + fread(fields->vy, sizeof(float), fields->yface_len, f) != + fields->yface_len) { + fclose(f); + fa_die("cannot read raw field payload"); + } + memcpy(fields->stage[0], fields->base, fields->scalar_len * sizeof(float)); + memcpy(fields->stage[1], fields->base, fields->scalar_len * sizeof(float)); + } else if (strcmp(mode, "analytic") == 0) { + fa_generate_analytic(cfg, fields); + } else { + fclose(f); + fa_die("unknown input mode"); + } + fclose(f); +} + +static void fa_write_csv(const char *path, const fa_config_t *cfg, + const fa_fields_t *fields) +{ + FILE *f = fopen(path, "wb"); + if (!f) { + fa_dief("cannot open output file: %s", path); + } + fprintf(f, "x,y,scalar\n"); + for (size_t j = 0; j < (size_t)cfg->ny + 2; ++j) { + const char *y_text = fa_center_coord_text(j, cfg->ny); + const float y = y_text ? strtof(y_text, NULL) : + fa_center_coord(j, cfg->ny, -FA_PI, FA_PI); + for (size_t i = 0; i < (size_t)cfg->nx + 2; ++i) { + const char *x_text = fa_center_coord_text(i, cfg->nx); + const float x = x_text ? strtof(x_text, NULL) : + fa_center_coord(i, cfg->nx, -FA_PI, FA_PI); + const size_t idx = fa_scalar_index(i, j, cfg->nx); + if (x_text && y_text) { + fprintf(f, "%s,%s,%.13g\n", x_text, y_text, + (double)fields->base[idx]); + } else { + fprintf(f, "%.12f,%.12f,%.13g\n", (double)x, (double)y, + (double)fields->base[idx]); + } + } + if (j + 1 < (size_t)cfg->ny + 2) { + fputc('\n', f); + } + } + fclose(f); +} + +static void fa_write_binary(const char *path, const fa_fields_t *fields) +{ + FILE *f = fopen(path, "wb"); + if (!f) { + fa_dief("cannot open output file: %s", path); + } + if (fwrite(fields->base, sizeof(float), fields->scalar_len, f) != + fields->scalar_len) { + fclose(f); + fa_die("cannot write binary field payload"); + } + fclose(f); +} + +static void fa_check_finite(const float *data, size_t n, const char *name) +{ + for (size_t i = 0; i < n; ++i) { + if (!isfinite(data[i])) { + fa_dief("%s contains non-finite values", name); + } + } +} + +static inline float fa_grad_face_line(const float *line, size_t stride, size_t n, size_t i, + float inv_h) +{ + if (i == 0) { + return (((-352.0f / 105.0f) * inv_h) * line[0 * stride] + + ((35.0f / 8.0f) * inv_h) * line[1 * stride] + + ((-35.0f / 24.0f) * inv_h) * line[2 * stride] + + ((21.0f / 40.0f) * inv_h) * line[3 * stride] + + ((-5.0f / 56.0f) * inv_h) * line[4 * stride]); + } + if (i == 1) { + return (((16.0f / 105.0f) * inv_h) * line[0 * stride] + + ((-31.0f / 24.0f) * inv_h) * line[1 * stride] + + ((29.0f / 24.0f) * inv_h) * line[2 * stride] + + ((-3.0f / 40.0f) * inv_h) * line[3 * stride] + + ((1.0f / 168.0f) * inv_h) * line[4 * stride]); + } + if (i + 3 < n) { + return (((1.0f / 24.0f) * inv_h) * line[(i - 1) * stride] + + ((-9.0f / 8.0f) * inv_h) * line[i * stride] + + ((9.0f / 8.0f) * inv_h) * line[(i + 1) * stride] + + ((-1.0f / 24.0f) * inv_h) * line[(i + 2) * stride]); + } + if (i == n - 3) { + return (((-1.0f / 168.0f) * inv_h) * line[(n - 5) * stride] + + ((3.0f / 40.0f) * inv_h) * line[(n - 4) * stride] + + ((-29.0f / 24.0f) * inv_h) * line[(n - 3) * stride] + + ((31.0f / 24.0f) * inv_h) * line[(n - 2) * stride] + + ((-16.0f / 105.0f) * inv_h) * line[(n - 1) * stride]); + } + return (((5.0f / 56.0f) * inv_h) * line[(n - 5) * stride] + + ((-21.0f / 40.0f) * inv_h) * line[(n - 4) * stride] + + ((35.0f / 24.0f) * inv_h) * line[(n - 3) * stride] + + ((-35.0f / 8.0f) * inv_h) * line[(n - 2) * stride] + + ((352.0f / 105.0f) * inv_h) * line[(n - 1) * stride]); +} + +static inline float fa_interp_c2f_line(const float *line, size_t stride, size_t n, size_t i) +{ + if (i == 0) { + return line[0 * stride]; + } + if (i == 1) { + return (-16.0f / 112.0f) * line[0 * stride] + + (70.0f / 112.0f) * line[1 * stride] + + (70.0f / 112.0f) * line[2 * stride] + + (-14.0f / 112.0f) * line[3 * stride] + + (2.0f / 112.0f) * line[4 * stride]; + } + if (i + 3 < n) { + return (-7.0f / 112.0f) * line[(i - 1) * stride] + + (63.0f / 112.0f) * line[i * stride] + + (63.0f / 112.0f) * line[(i + 1) * stride] + + (-7.0f / 112.0f) * line[(i + 2) * stride]; + } + if (i == n - 3) { + return (2.0f / 112.0f) * line[(n - 5) * stride] + + (-14.0f / 112.0f) * line[(n - 4) * stride] + + (70.0f / 112.0f) * line[(n - 3) * stride] + + (70.0f / 112.0f) * line[(n - 2) * stride] + + (-16.0f / 112.0f) * line[(n - 1) * stride]; + } + if (i == n - 2) { + return line[(n - 1) * stride]; + } + fa_die("invalid center-to-face stencil index"); + return 0.0f; +} + +static inline float fa_div_face_line(const float *line, size_t stride, size_t n, size_t i, + float inv_h) +{ + if (i == 0) { + return (((-11.0f / 12.0f) * inv_h) * line[0 * stride] + + ((17.0f / 24.0f) * inv_h) * line[1 * stride] + + ((3.0f / 8.0f) * inv_h) * line[2 * stride] + + ((-5.0f / 24.0f) * inv_h) * line[3 * stride] + + ((1.0f / 24.0f) * inv_h) * line[4 * stride]); + } + if (i + 2 < n) { + return (((1.0f / 24.0f) * inv_h) * line[(i - 1) * stride] + + ((-9.0f / 8.0f) * inv_h) * line[i * stride] + + ((9.0f / 8.0f) * inv_h) * line[(i + 1) * stride] + + ((-1.0f / 24.0f) * inv_h) * line[(i + 2) * stride]); + } + return (((-1.0f / 24.0f) * inv_h) * line[(n - 5) * stride] + + ((5.0f / 24.0f) * inv_h) * line[(n - 4) * stride] + + ((-3.0f / 8.0f) * inv_h) * line[(n - 3) * stride] + + ((-17.0f / 24.0f) * inv_h) * line[(n - 2) * stride] + + ((11.0f / 12.0f) * inv_h) * line[(n - 1) * stride]); +} + +static void fa_rhs_cpu(const fa_config_t *cfg, const fa_fields_t *fields, + const float *stage_in, float *rhs_out) +{ + const float inv_dx = (float)cfg->nx / (2.0f * FA_PI); + const float inv_dy = (float)cfg->ny / (2.0f * FA_PI); + const size_t sx = (size_t)cfg->nx + 2; + const size_t sy = (size_t)cfg->ny + 2; + + float *diff_x = calloc(fields->xface_len, sizeof(float)); + float *adv_x = calloc(fields->xface_len, sizeof(float)); + float *diff_y = calloc(fields->yface_len, sizeof(float)); + float *adv_y = calloc(fields->yface_len, sizeof(float)); + if (!diff_x || !adv_x || !diff_y || !adv_y) { + fa_die("cannot allocate RHS scratch"); + } + + for (size_t j = 0; j < sy; ++j) { + const float *line = stage_in + j * sx; + for (size_t i = 0; i < sx - 1; ++i) { + const size_t idx = fa_xface_index(i, j, cfg->nx); + diff_x[idx] = FA_DIFFUSION * fa_grad_face_line(line, 1, sx, i, inv_dx); + adv_x[idx] = fields->vx[idx] * fa_interp_c2f_line(line, 1, sx, i); + } + } + + for (size_t i = 0; i < sx; ++i) { + const float *line = stage_in + i; + for (size_t j = 0; j < sy - 1; ++j) { + const size_t idx = fa_yface_index(i, j, cfg->nx); + diff_y[idx] = FA_DIFFUSION * fa_grad_face_line(line, sx, sy, j, inv_dy); + adv_y[idx] = fields->vy[idx] * fa_interp_c2f_line(line, sx, sy, j); + } + } + + memset(rhs_out, 0, fields->scalar_len * sizeof(float)); + for (size_t j = 0; j < (size_t)cfg->ny; ++j) { + for (size_t i = 0; i < (size_t)cfg->nx; ++i) { + const float diff = fa_div_face_line(diff_x + (j + 1) * ((size_t)cfg->nx + 1), + 1, (size_t)cfg->nx + 1, i, inv_dx) + + fa_div_face_line(diff_y + (i + 1), sx, (size_t)cfg->ny + 1, j, inv_dy); + const float adv = fa_div_face_line(adv_x + (j + 1) * ((size_t)cfg->nx + 1), + 1, (size_t)cfg->nx + 1, i, inv_dx) + + fa_div_face_line(adv_y + (i + 1), sx, (size_t)cfg->ny + 1, j, inv_dy); + rhs_out[fa_scalar_index(i + 1, j + 1, cfg->nx)] = diff - adv; + } + } + + free(diff_x); + free(adv_x); + free(diff_y); + free(adv_y); +} + +static void fa_rk4_step_cpu(const fa_config_t *cfg, fa_fields_t *fields) +{ + const float dt = cfg->dt; + const float dt_half = 0.5f * dt; + const float dt_over_6 = dt / 6.0f; + float *rhs = calloc(fields->scalar_len, sizeof(float)); + if (!rhs) { + fa_die("cannot allocate RHS buffer"); + } + + for (uint32_t stage = 1; stage <= 4; ++stage) { + const float *stage_in = fields->stage[(stage - 1) & 1u]; + float *stage_out = fields->stage[stage & 1u]; + fa_rhs_cpu(cfg, fields, stage_in, rhs); + + for (size_t j = 0; j < (size_t)cfg->ny + 2; ++j) { + for (size_t i = 0; i < (size_t)cfg->nx + 2; ++i) { + const size_t idx = fa_scalar_index(i, j, cfg->nx); + if (i == 0 || j == 0 || i == (size_t)cfg->nx + 1 || + j == (size_t)cfg->ny + 1) { + stage_out[idx] = fields->base[idx]; + if (stage == 4) { + fields->base[idx] = fields->base[idx]; + } + continue; + } + if (stage == 1) { + fields->accum[idx] = rhs[idx]; + stage_out[idx] = fields->base[idx] + dt_half * rhs[idx]; + } else if (stage == 2) { + fields->accum[idx] = fields->accum[idx] + 2.0f * rhs[idx]; + stage_out[idx] = fields->base[idx] + dt_half * rhs[idx]; + } else if (stage == 3) { + fields->accum[idx] = fields->accum[idx] + 2.0f * rhs[idx]; + stage_out[idx] = fields->base[idx] + dt * rhs[idx]; + } else { + fields->accum[idx] = fields->accum[idx] + rhs[idx]; + fields->base[idx] = fields->base[idx] + dt_over_6 * fields->accum[idx]; + stage_out[idx] = fields->base[idx]; + } + } + } + } + + free(rhs); +} + +static id fa_create_device(void) +{ + id device = MTLCreateSystemDefaultDevice(); + if (!device) { + fa_die("no Metal device available"); + } + return device; +} + +static NSString *fa_read_text_file(NSString *path) +{ + NSError *error = nil; + NSString *content = [NSString stringWithContentsOfFile:path + encoding:NSUTF8StringEncoding + error:&error]; + if (!content) { + fprintf(stderr, "cannot read %s: %s\n", path.UTF8String, + error.localizedDescription.UTF8String); + exit(1); + } + return content; +} + +typedef struct { + uint32_t nx; + uint32_t ny; + uint32_t stage; + uint32_t unused0; + float inv_dx; + float inv_dy; + float alpha; + float dt_over_6; + float weight; + float unused1; + float unused2; + float unused3; +} fa_metal_params_t; + +_Static_assert(sizeof(fa_metal_params_t) == 48, "fa_metal_params_t must be 48 bytes"); +_Static_assert(offsetof(fa_metal_params_t, nx) == 0, "fa_metal_params_t.nx"); +_Static_assert(offsetof(fa_metal_params_t, ny) == 4, "fa_metal_params_t.ny"); +_Static_assert(offsetof(fa_metal_params_t, stage) == 8, "fa_metal_params_t.stage"); +_Static_assert(offsetof(fa_metal_params_t, inv_dx) == 16, "fa_metal_params_t.inv_dx"); +_Static_assert(offsetof(fa_metal_params_t, inv_dy) == 20, "fa_metal_params_t.inv_dy"); +_Static_assert(offsetof(fa_metal_params_t, alpha) == 24, "fa_metal_params_t.alpha"); +_Static_assert(offsetof(fa_metal_params_t, dt_over_6) == 28, "fa_metal_params_t.dt_over_6"); +_Static_assert(offsetof(fa_metal_params_t, weight) == 32, "fa_metal_params_t.weight"); + +static NSString *fa_metal_source_path(void) +{ + NSString *cwd = [[NSFileManager defaultManager] currentDirectoryPath]; + return [cwd stringByAppendingPathComponent:@"formal_advection.metal"]; +} + +static bool fa_metal_timestamps_valid(const id cb) +{ + return cb.GPUStartTime > 0.0 && cb.GPUEndTime > cb.GPUStartTime; +} + +static void fa_run_metal(const fa_options_t *opt, fa_fields_t *fields, + fa_run_stats_t *stats) +{ + stats->jit_seconds = 0.0; + stats->solve_seconds = 0.0; + stats->gpu_busy = 0.0; + stats->final_all_finite = false; + + double jit_start = fa_now_seconds(); + id device = fa_create_device(); + NSError *error = nil; + NSString *source = fa_read_text_file(fa_metal_source_path()); + MTLCompileOptions *compile_options = [MTLCompileOptions new]; + id library = [device newLibraryWithSource:source + options:compile_options + error:&error]; + if (!library) { + fprintf(stderr, "Metal compile failed: %s\n", + error.localizedDescription.UTF8String); + exit(1); + } + + const char *flux_name = opt->index_mode == FA_INDEX_LINEAR ? + "compute_fluxes_linear" : "compute_fluxes_2d"; + const char *update_name = opt->index_mode == FA_INDEX_LINEAR ? + "update_stage_linear" : "update_stage_2d"; + id flux_fn = [library newFunctionWithName:[NSString stringWithUTF8String:flux_name]]; + id update_fn = [library newFunctionWithName:[NSString stringWithUTF8String:update_name]]; + if (!flux_fn || !update_fn) { + fa_die("cannot find Metal kernels"); + } + id flux_pso = + [device newComputePipelineStateWithFunction:flux_fn error:&error]; + id update_pso = + [device newComputePipelineStateWithFunction:update_fn error:&error]; + if (!flux_pso || !update_pso) { + fprintf(stderr, "cannot create Metal pipelines: %s\n", + error.localizedDescription.UTF8String); + exit(1); + } + stats->jit_seconds = fa_now_seconds() - jit_start; + + id queue = [device newCommandQueue]; + if (!queue) { + fa_die("cannot create Metal command queue"); + } + + if (opt->tg_supplied) { + if (opt->tg_x == 0 || opt->tg_y == 0) { + fa_die("threadgroup dimensions must be positive"); + } + if (opt->index_mode == FA_INDEX_LINEAR && opt->tg_y != 1) { + fa_die("linear index mode requires 1D threadgroups"); + } + if ((uint64_t)opt->tg_x * (uint64_t)opt->tg_y > + (uint64_t)flux_pso.maxTotalThreadsPerThreadgroup) { + fa_die("threadgroup exceeds pipeline limit"); + } + } + + float *base_backup = malloc(fields->scalar_len * sizeof(float)); + float *stage0_backup = malloc(fields->scalar_len * sizeof(float)); + float *stage1_backup = malloc(fields->scalar_len * sizeof(float)); + float *accum_backup = malloc(fields->scalar_len * sizeof(float)); + if (!base_backup || !stage0_backup || !stage1_backup || !accum_backup) { + fa_die("cannot allocate warmup backups"); + } + memcpy(base_backup, fields->base, fields->scalar_len * sizeof(float)); + memcpy(stage0_backup, fields->stage[0], fields->scalar_len * sizeof(float)); + memcpy(stage1_backup, fields->stage[1], fields->scalar_len * sizeof(float)); + memcpy(accum_backup, fields->accum, fields->scalar_len * sizeof(float)); + + id base = [device newBufferWithBytes:fields->base + length:fields->scalar_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id stage0 = [device newBufferWithBytes:fields->stage[0] + length:fields->scalar_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id stage1 = [device newBufferWithBytes:fields->stage[1] + length:fields->scalar_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id accum = [device newBufferWithBytes:fields->accum + length:fields->scalar_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id vx = [device newBufferWithBytes:fields->vx + length:fields->xface_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id vy = [device newBufferWithBytes:fields->vy + length:fields->yface_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id diff_x = [device newBufferWithLength:fields->xface_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id adv_x = [device newBufferWithLength:fields->xface_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id diff_y = [device newBufferWithLength:fields->yface_len * sizeof(float) + options:MTLResourceStorageModeShared]; + id adv_y = [device newBufferWithLength:fields->yface_len * sizeof(float) + options:MTLResourceStorageModeShared]; + if (!base || !stage0 || !stage1 || !accum || !vx || !vy || + !diff_x || !adv_x || !diff_y || !adv_y) { + fa_die("cannot allocate Metal buffers"); + } + + const MTLSize tg = MTLSizeMake(opt->tg_x, opt->tg_y, 1); + const MTLSize grid = opt->index_mode == FA_INDEX_LINEAR ? + MTLSizeMake((NSUInteger)(opt->cfg.nx + 2) * (NSUInteger)(opt->cfg.ny + 2), 1, 1) : + MTLSizeMake((NSUInteger)opt->cfg.nx + 2, + (NSUInteger)opt->cfg.ny + 2, 1); + const float dt_half = 0.5f * opt->cfg.dt; + const float dt_over_6 = opt->cfg.dt / 6.0f; + const uint32_t steps_per_cb = 8; + id in_buf = stage0; + id out_buf = stage1; + + { + id warmup_cb = [queue commandBuffer]; + id warmup_enc = [warmup_cb computeCommandEncoder]; + for (uint32_t stage = 1; stage <= 4; ++stage) { + fa_metal_params_t warmup_params = { + .nx = opt->cfg.nx, + .ny = opt->cfg.ny, + .stage = stage, + .inv_dx = (float)opt->cfg.nx / (2.0f * FA_PI), + .inv_dy = (float)opt->cfg.ny / (2.0f * FA_PI), + .alpha = (stage == 1 || stage == 2) ? dt_half : + (stage == 3 ? opt->cfg.dt : 0.0f), + .dt_over_6 = dt_over_6, + .weight = (stage == 1 || stage == 4) ? 1.0f : 2.0f, + }; + [warmup_enc setComputePipelineState:flux_pso]; + [warmup_enc setBuffer:in_buf offset:0 atIndex:0]; + [warmup_enc setBuffer:vx offset:0 atIndex:1]; + [warmup_enc setBuffer:vy offset:0 atIndex:2]; + [warmup_enc setBuffer:diff_x offset:0 atIndex:3]; + [warmup_enc setBuffer:adv_x offset:0 atIndex:4]; + [warmup_enc setBuffer:diff_y offset:0 atIndex:5]; + [warmup_enc setBuffer:adv_y offset:0 atIndex:6]; + [warmup_enc setBytes:&warmup_params length:sizeof(warmup_params) atIndex:7]; + [warmup_enc dispatchThreads:grid threadsPerThreadgroup:tg]; + + [warmup_enc setComputePipelineState:update_pso]; + [warmup_enc setBuffer:base offset:0 atIndex:0]; + [warmup_enc setBuffer:in_buf offset:0 atIndex:1]; + [warmup_enc setBuffer:out_buf offset:0 atIndex:2]; + [warmup_enc setBuffer:accum offset:0 atIndex:3]; + [warmup_enc setBuffer:diff_x offset:0 atIndex:4]; + [warmup_enc setBuffer:adv_x offset:0 atIndex:5]; + [warmup_enc setBuffer:diff_y offset:0 atIndex:6]; + [warmup_enc setBuffer:adv_y offset:0 atIndex:7]; + [warmup_enc setBytes:&warmup_params length:sizeof(warmup_params) atIndex:8]; + [warmup_enc dispatchThreads:grid threadsPerThreadgroup:tg]; + + id tmp = in_buf; + in_buf = out_buf; + out_buf = tmp; + } + [warmup_enc endEncoding]; + [warmup_cb commit]; + [warmup_cb waitUntilCompleted]; + if (warmup_cb.error) { + fprintf(stderr, "Metal warmup failed: %s\n", + warmup_cb.error.localizedDescription.UTF8String); + exit(1); + } + if (!fa_metal_timestamps_valid(warmup_cb)) { + fa_die("Metal warmup produced invalid timestamps"); + } + + memcpy(base.contents, base_backup, fields->scalar_len * sizeof(float)); + memcpy(stage0.contents, stage0_backup, fields->scalar_len * sizeof(float)); + memcpy(stage1.contents, stage1_backup, fields->scalar_len * sizeof(float)); + memcpy(accum.contents, accum_backup, fields->scalar_len * sizeof(float)); + in_buf = stage0; + out_buf = stage1; + } + + double solve_start = fa_now_seconds(); + __strong id pending[2] = { nil, nil }; + uint32_t pending_count = 0; + double gpu_busy = 0.0; + double gpu_cursor = 0.0; + bool have_gpu_cursor = false; + for (uint32_t step_base = 0; step_base < opt->cfg.steps; ) { + const uint32_t batch_steps = (opt->cfg.steps - step_base < steps_per_cb) ? + (opt->cfg.steps - step_base) : steps_per_cb; + @autoreleasepool { + id cb = [queue commandBuffer]; + id enc = [cb computeCommandEncoder]; + for (uint32_t local_step = 0; local_step < batch_steps; ++local_step) { + for (uint32_t stage = 1; stage <= 4; ++stage) { + fa_metal_params_t params = { + .nx = opt->cfg.nx, + .ny = opt->cfg.ny, + .stage = stage, + .inv_dx = (float)opt->cfg.nx / (2.0f * FA_PI), + .inv_dy = (float)opt->cfg.ny / (2.0f * FA_PI), + .alpha = (stage == 1 || stage == 2) ? dt_half : + (stage == 3 ? opt->cfg.dt : 0.0f), + .dt_over_6 = dt_over_6, + .weight = (stage == 1 || stage == 4) ? 1.0f : 2.0f, + }; + [enc setComputePipelineState:flux_pso]; + [enc setBuffer:in_buf offset:0 atIndex:0]; + [enc setBuffer:vx offset:0 atIndex:1]; + [enc setBuffer:vy offset:0 atIndex:2]; + [enc setBuffer:diff_x offset:0 atIndex:3]; + [enc setBuffer:adv_x offset:0 atIndex:4]; + [enc setBuffer:diff_y offset:0 atIndex:5]; + [enc setBuffer:adv_y offset:0 atIndex:6]; + [enc setBytes:¶ms length:sizeof(params) atIndex:7]; + [enc dispatchThreads:grid threadsPerThreadgroup:tg]; + + [enc setComputePipelineState:update_pso]; + [enc setBuffer:base offset:0 atIndex:0]; + [enc setBuffer:in_buf offset:0 atIndex:1]; + [enc setBuffer:out_buf offset:0 atIndex:2]; + [enc setBuffer:accum offset:0 atIndex:3]; + [enc setBuffer:diff_x offset:0 atIndex:4]; + [enc setBuffer:adv_x offset:0 atIndex:5]; + [enc setBuffer:diff_y offset:0 atIndex:6]; + [enc setBuffer:adv_y offset:0 atIndex:7]; + [enc setBytes:¶ms length:sizeof(params) atIndex:8]; + [enc dispatchThreads:grid threadsPerThreadgroup:tg]; + + id tmp = in_buf; + in_buf = out_buf; + out_buf = tmp; + } + } + [enc endEncoding]; + [cb commit]; + pending[pending_count++] = cb; + if (pending_count == 2) { + id oldest = pending[0]; + [oldest waitUntilCompleted]; + if (oldest.error) { + fprintf(stderr, "Metal solve failed: %s\n", + oldest.error.localizedDescription.UTF8String); + exit(1); + } + if (!fa_metal_timestamps_valid(oldest)) { + fa_die("Metal solve produced invalid timestamps"); + } + if (!have_gpu_cursor) { + gpu_busy += oldest.GPUEndTime - oldest.GPUStartTime; + gpu_cursor = oldest.GPUEndTime; + have_gpu_cursor = true; + } else if (oldest.GPUStartTime > gpu_cursor) { + gpu_busy += oldest.GPUEndTime - oldest.GPUStartTime; + gpu_cursor = oldest.GPUEndTime; + } else if (oldest.GPUEndTime > gpu_cursor) { + gpu_busy += oldest.GPUEndTime - gpu_cursor; + gpu_cursor = oldest.GPUEndTime; + } + pending[0] = pending[1]; + pending[1] = nil; + pending_count = 1; + } + } + step_base += batch_steps; + } + + for (uint32_t i = 0; i < pending_count; ++i) { + id cb = pending[i]; + [cb waitUntilCompleted]; + if (cb.error) { + fprintf(stderr, "Metal solve failed: %s\n", + cb.error.localizedDescription.UTF8String); + exit(1); + } + if (!fa_metal_timestamps_valid(cb)) { + fa_die("Metal solve produced invalid timestamps"); + } + if (!have_gpu_cursor) { + gpu_busy += cb.GPUEndTime - cb.GPUStartTime; + gpu_cursor = cb.GPUEndTime; + have_gpu_cursor = true; + } else if (cb.GPUStartTime > gpu_cursor) { + gpu_busy += cb.GPUEndTime - cb.GPUStartTime; + gpu_cursor = cb.GPUEndTime; + } else if (cb.GPUEndTime > gpu_cursor) { + gpu_busy += cb.GPUEndTime - gpu_cursor; + gpu_cursor = cb.GPUEndTime; + } + } + memcpy(fields->base, base.contents, fields->scalar_len * sizeof(float)); + stats->solve_seconds = fa_now_seconds() - solve_start; + stats->gpu_busy = gpu_busy; + fa_check_finite(fields->base, fields->scalar_len, "final scalar"); + stats->final_all_finite = true; + + free(base_backup); + free(stage0_backup); + free(stage1_backup); + free(accum_backup); +} + +static void fa_run_cpu(const fa_options_t *opt, fa_fields_t *fields, + fa_run_stats_t *stats) +{ + stats->jit_seconds = 0.0; + stats->gpu_busy = 0.0; + double t0 = fa_now_seconds(); + for (uint32_t step = 0; step < opt->cfg.steps; ++step) { + fa_rk4_step_cpu(&opt->cfg, fields); + } + stats->solve_seconds = fa_now_seconds() - t0; + fa_check_finite(fields->base, fields->scalar_len, "final scalar"); + stats->final_all_finite = true; +} + +static void fa_usage(const char *argv0) +{ + fprintf(stderr, + "usage: %s [--mode cpu|metal] [--cells nx,ny] [--steps n] [--dt dt]\n" + " [--input path] [--output path] [--dump [path]] [--dump-bin [path]]\n" + " [--index-mode 2d|linear] [--tg x[,y]]\n", + argv0); +} + +static bool fa_parse_cells(const char *raw, fa_config_t *cfg) +{ + const char *comma = strchr(raw, ','); + if (!comma) { + return false; + } + size_t left_len = (size_t)(comma - raw); + size_t right_len = strlen(comma + 1); + if (left_len == 0 || right_len == 0) { + return false; + } + char left[32]; + char right[32]; + if (left_len >= sizeof(left) || right_len >= sizeof(right)) { + return false; + } + memcpy(left, raw, left_len); + left[left_len] = '\0'; + memcpy(right, comma + 1, right_len + 1); + cfg->nx = fa_parse_u32_strict(left, "--cells"); + cfg->ny = fa_parse_u32_strict(right, "--cells"); + return true; +} + +static bool fa_parse_threadgroup(const char *raw, uint32_t *x, uint32_t *y) +{ + const char *comma = strchr(raw, ','); + if (!comma) { + *x = fa_parse_u32_strict(raw, "--tg"); + *y = 1; + return true; + } + size_t left_len = (size_t)(comma - raw); + size_t right_len = strlen(comma + 1); + if (left_len == 0 || right_len == 0 || left_len >= 32 || right_len >= 32) { + return false; + } + char left[32]; + char right[32]; + memcpy(left, raw, left_len); + left[left_len] = '\0'; + memcpy(right, comma + 1, right_len + 1); + *x = fa_parse_u32_strict(left, "--tg"); + *y = fa_parse_u32_strict(right, "--tg"); + return true; +} + +int main(int argc, char **argv) +{ + @autoreleasepool { + fa_options_t opt = { + .metal = false, + .dump_csv = false, + .dump_bin = false, + .input_path = NULL, + .output_path = NULL, + .binary_output_path = NULL, + .cfg = { + .nx = 51, + .ny = 51, + .steps = 500, + .dt = 1.0e-4f, + }, + .dt_supplied = false, + .tg_x = 32, + .tg_y = 4, + .tg_supplied = false, + .index_mode_supplied = false, + .index_mode = FA_INDEX_DIRECT_2D, + }; + + for (int i = 1; i < argc; ++i) { + if (strcmp(argv[i], "--help") == 0 || strcmp(argv[i], "-h") == 0) { + fa_usage(argv[0]); + return 0; + } else if (strcmp(argv[i], "--mode") == 0 && i + 1 < argc) { + ++i; + if (strcmp(argv[i], "cpu") == 0) { + opt.metal = false; + } else if (strcmp(argv[i], "metal") == 0) { + opt.metal = true; + } else { + fa_die("mode must be cpu or metal"); + } + } else if (strcmp(argv[i], "--cells") == 0 && i + 1 < argc) { + if (!fa_parse_cells(argv[++i], &opt.cfg)) { + fa_die("cannot parse --cells nx,ny"); + } + } else if (strcmp(argv[i], "--steps") == 0 && i + 1 < argc) { + opt.cfg.steps = fa_parse_u32_strict(argv[++i], "--steps"); + } else if (strcmp(argv[i], "--dt") == 0 && i + 1 < argc) { + opt.cfg.dt = fa_parse_float_strict(argv[++i], "--dt"); + opt.dt_supplied = true; + } else if (strcmp(argv[i], "--input") == 0 && i + 1 < argc) { + opt.input_path = argv[++i]; + } else if (strcmp(argv[i], "--output") == 0 && i + 1 < argc) { + opt.output_path = argv[++i]; + opt.dump_csv = true; + } else if (strcmp(argv[i], "--dump-bin") == 0) { + opt.dump_bin = true; + if (i + 1 < argc && argv[i + 1][0] != '-') { + opt.binary_output_path = argv[++i]; + } + } else if (strcmp(argv[i], "--index-mode") == 0 && i + 1 < argc) { + ++i; + if (strcmp(argv[i], "2d") == 0) { + opt.index_mode = FA_INDEX_DIRECT_2D; + } else if (strcmp(argv[i], "linear") == 0) { + opt.index_mode = FA_INDEX_LINEAR; + } else { + fa_die("index mode must be 2d or linear"); + } + opt.index_mode_supplied = true; + } else if (strcmp(argv[i], "--tg") == 0 && i + 1 < argc) { + if (!fa_parse_threadgroup(argv[++i], &opt.tg_x, &opt.tg_y)) { + fa_die("cannot parse --tg x[,y]"); + } + opt.tg_supplied = true; + } else if (strcmp(argv[i], "--dump") == 0) { + opt.dump_csv = true; + if (i + 1 < argc && argv[i + 1][0] != '-') { + opt.output_path = argv[++i]; + } + } else { + fa_dief("unknown argument: %s", argv[i]); + } + } + + if (!opt.index_mode_supplied) { + const char *env_index_mode = getenv("FA_INDEX_MODE"); + if (env_index_mode) { + if (strcmp(env_index_mode, "2d") == 0) { + opt.index_mode = FA_INDEX_DIRECT_2D; + } else if (strcmp(env_index_mode, "linear") == 0) { + opt.index_mode = FA_INDEX_LINEAR; + } else { + fa_die("FA_INDEX_MODE must be 2d or linear"); + } + } + } + if (!opt.tg_supplied) { + const char *env_tg = getenv("FA_TG"); + if (env_tg) { + if (!fa_parse_threadgroup(env_tg, &opt.tg_x, &opt.tg_y)) { + fa_die("FA_TG must be x[,y]"); + } + opt.tg_supplied = true; + } + } + if (opt.index_mode == FA_INDEX_LINEAR && !opt.tg_supplied) { + opt.tg_x = 64; + opt.tg_y = 1; + } + fa_validate_cfg(&opt.cfg, opt.dt_supplied); + fa_validate_working_set(&opt.cfg); + + const double t_process_start = fa_now_seconds(); + double t_setup_start = t_process_start; + fa_fields_t fields; + if (opt.input_path) { + fa_read_raw_input(opt.input_path, &opt.cfg, &fields); + } else { + fa_fields_alloc(&fields, &opt.cfg); + fa_generate_analytic(&opt.cfg, &fields); + } + fa_check_finite(fields.base, fields.scalar_len, "initial scalar"); + fa_check_finite(fields.vx, fields.xface_len, "initial vx"); + fa_check_finite(fields.vy, fields.yface_len, "initial vy"); + const double t_setup_end = fa_now_seconds(); + fa_run_stats_t stats; + if (opt.metal) { + fa_run_metal(&opt, &fields, &stats); + } else { + fa_run_cpu(&opt, &fields, &stats); + } + + if (opt.dump_csv) { + if (!opt.output_path) { + opt.output_path = "final.csv"; + } + fa_write_csv(opt.output_path, &opt.cfg, &fields); + } + if (opt.dump_bin) { + if (!opt.binary_output_path) { + opt.binary_output_path = "final.bin"; + } + fa_write_binary(opt.binary_output_path, &fields); + } + if (!stats.final_all_finite) { + fa_die("final field contains non-finite values"); + } + const double t_process_end = fa_now_seconds(); + fprintf(stderr, + "[startup] setup=%.6f jit=%.6f solve=%.6f total_elapsed=%.6f index_mode=%s tg=%u,%u\n", + t_setup_end - t_setup_start, + stats.jit_seconds, + stats.solve_seconds, + t_process_end - t_process_start, + opt.index_mode == FA_INDEX_LINEAR ? "linear" : "2d", + opt.tg_x, opt.tg_y); + fprintf(stdout, + "RESULT mode=%s solve_wall=%.6f gpu_busy=%.6f cells=%u,%u steps=%u dt=%.9g index_mode=%s tg=%u,%u\n", + opt.metal ? "metal" : "cpu", + stats.solve_seconds, + stats.gpu_busy, + opt.cfg.nx, opt.cfg.ny, opt.cfg.steps, opt.cfg.dt, + opt.index_mode == FA_INDEX_LINEAR ? "linear" : "2d", + opt.tg_x, opt.tg_y); + fa_fields_free(&fields); + return 0; + } +} diff --git a/benchmarks/2d-advection-metal/formal_advection.metal b/benchmarks/2d-advection-metal/formal_advection.metal new file mode 100644 index 0000000..0e49203 --- /dev/null +++ b/benchmarks/2d-advection-metal/formal_advection.metal @@ -0,0 +1,254 @@ +#include +using namespace metal; + +struct fa_params_t { + uint nx; + uint ny; + uint stage; + uint unused0; + float inv_dx; + float inv_dy; + float alpha; + float dt_over_6; + float weight; + float unused1; + float unused2; + float unused3; +}; + +inline float fa_grad_face_line(const device float *line, uint stride, uint n, uint i, + float inv_h) +{ + if (i == 0) { + return (((-352.0f / 105.0f) * inv_h) * line[0 * stride] + + ((35.0f / 8.0f) * inv_h) * line[1 * stride] + + ((-35.0f / 24.0f) * inv_h) * line[2 * stride] + + ((21.0f / 40.0f) * inv_h) * line[3 * stride] + + ((-5.0f / 56.0f) * inv_h) * line[4 * stride]); + } + if (i == 1) { + return (((16.0f / 105.0f) * inv_h) * line[0 * stride] + + ((-31.0f / 24.0f) * inv_h) * line[1 * stride] + + ((29.0f / 24.0f) * inv_h) * line[2 * stride] + + ((-3.0f / 40.0f) * inv_h) * line[3 * stride] + + ((1.0f / 168.0f) * inv_h) * line[4 * stride]); + } + if (i + 3 < n) { + return (((1.0f / 24.0f) * inv_h) * line[(i - 1) * stride] + + ((-9.0f / 8.0f) * inv_h) * line[i * stride] + + ((9.0f / 8.0f) * inv_h) * line[(i + 1) * stride] + + ((-1.0f / 24.0f) * inv_h) * line[(i + 2) * stride]); + } + if (i == n - 3) { + return (((-1.0f / 168.0f) * inv_h) * line[(n - 5) * stride] + + ((3.0f / 40.0f) * inv_h) * line[(n - 4) * stride] + + ((-29.0f / 24.0f) * inv_h) * line[(n - 3) * stride] + + ((31.0f / 24.0f) * inv_h) * line[(n - 2) * stride] + + ((-16.0f / 105.0f) * inv_h) * line[(n - 1) * stride]); + } + return (((5.0f / 56.0f) * inv_h) * line[(n - 5) * stride] + + ((-21.0f / 40.0f) * inv_h) * line[(n - 4) * stride] + + ((35.0f / 24.0f) * inv_h) * line[(n - 3) * stride] + + ((-35.0f / 8.0f) * inv_h) * line[(n - 2) * stride] + + ((352.0f / 105.0f) * inv_h) * line[(n - 1) * stride]); +} + +inline float fa_interp_c2f_line(const device float *line, uint stride, uint n, uint i) +{ + if (i == 0) { + return line[0 * stride]; + } + if (i == 1) { + return (-16.0f / 112.0f) * line[0 * stride] + + (70.0f / 112.0f) * line[1 * stride] + + (70.0f / 112.0f) * line[2 * stride] + + (-14.0f / 112.0f) * line[3 * stride] + + (2.0f / 112.0f) * line[4 * stride]; + } + if (i + 3 < n) { + return (-7.0f / 112.0f) * line[(i - 1) * stride] + + (63.0f / 112.0f) * line[i * stride] + + (63.0f / 112.0f) * line[(i + 1) * stride] + + (-7.0f / 112.0f) * line[(i + 2) * stride]; + } + if (i == n - 3) { + return (2.0f / 112.0f) * line[(n - 5) * stride] + + (-14.0f / 112.0f) * line[(n - 4) * stride] + + (70.0f / 112.0f) * line[(n - 3) * stride] + + (70.0f / 112.0f) * line[(n - 2) * stride] + + (-16.0f / 112.0f) * line[(n - 1) * stride]; + } + return line[(n - 1) * stride]; +} + +inline float fa_div_face_line(const device float *line, uint stride, uint n, uint i, + float inv_h) +{ + if (i == 0) { + return (((-11.0f / 12.0f) * inv_h) * line[0 * stride] + + ((17.0f / 24.0f) * inv_h) * line[1 * stride] + + ((3.0f / 8.0f) * inv_h) * line[2 * stride] + + ((-5.0f / 24.0f) * inv_h) * line[3 * stride] + + ((1.0f / 24.0f) * inv_h) * line[4 * stride]); + } + if (i + 2 < n) { + return (((1.0f / 24.0f) * inv_h) * line[(i - 1) * stride] + + ((-9.0f / 8.0f) * inv_h) * line[i * stride] + + ((9.0f / 8.0f) * inv_h) * line[(i + 1) * stride] + + ((-1.0f / 24.0f) * inv_h) * line[(i + 2) * stride]); + } + return (((-1.0f / 24.0f) * inv_h) * line[(n - 5) * stride] + + ((5.0f / 24.0f) * inv_h) * line[(n - 4) * stride] + + ((-3.0f / 8.0f) * inv_h) * line[(n - 3) * stride] + + ((-17.0f / 24.0f) * inv_h) * line[(n - 2) * stride] + + ((11.0f / 12.0f) * inv_h) * line[(n - 1) * stride]); +} + +inline uint fa_scalar_index(uint i, uint j, uint nx) +{ + return i + (nx + 2u) * j; +} + +inline uint fa_xface_index(uint i, uint j, uint nx) +{ + return i + (nx + 1u) * j; +} + +inline uint fa_yface_index(uint i, uint j, uint nx) +{ + return i + (nx + 2u) * j; +} + +inline void fa_compute_fluxes_body(device const float *stage, + device const float *vx, device const float *vy, device float *diff_x, + device float *adv_x, device float *diff_y, device float *adv_y, + constant fa_params_t &p, uint i, uint j) +{ + const uint nx = p.nx; + const uint ny = p.ny; + const uint sx = nx + 2u; + const uint sy = ny + 2u; + + if (i <= nx && j <= ny + 1u) { + const device float *line = stage + j * sx; + const uint idx = fa_xface_index(i, j, nx); + diff_x[idx] = 0.5f * fa_grad_face_line(line, 1u, sx, i, p.inv_dx); + adv_x[idx] = vx[idx] * fa_interp_c2f_line(line, 1u, sx, i); + } + + if (i <= nx + 1u && j <= ny) { + const device float *line = stage + i; + const uint idx = fa_yface_index(i, j, nx); + diff_y[idx] = 0.5f * fa_grad_face_line(line, sx, sy, j, p.inv_dy); + adv_y[idx] = vy[idx] * fa_interp_c2f_line(line, sx, sy, j); + } +} + +inline void fa_update_stage_body(device float *base, device const float *stage_in, + device float *stage_out, device float *accum, device const float *diff_x, + device const float *adv_x, device const float *diff_y, + device const float *adv_y, constant fa_params_t &p, uint i, uint j) +{ + const uint nx = p.nx; + const uint ny = p.ny; + const uint stage = p.stage; + const uint scalar_idx = fa_scalar_index(i, j, nx); + + if (i == 0u || j == 0u || i == nx + 1u || j == ny + 1u) { + stage_out[scalar_idx] = base[scalar_idx]; + return; + } + + const uint x_idx = i - 1u; + const uint y_idx = j - 1u; + const device float *x_line_diff = diff_x + j * (nx + 1u); + const device float *x_line_adv = adv_x + j * (nx + 1u); + const device float *y_line_diff = diff_y + i; + const device float *y_line_adv = adv_y + i; + + const float diffx = fa_div_face_line(x_line_diff, 1u, nx + 1u, x_idx, p.inv_dx); + const float diffy = fa_div_face_line(y_line_diff, nx + 2u, ny + 1u, y_idx, p.inv_dy); + const float advx = fa_div_face_line(x_line_adv, 1u, nx + 1u, x_idx, p.inv_dx); + const float advy = fa_div_face_line(y_line_adv, nx + 2u, ny + 1u, y_idx, p.inv_dy); + const float rhs = (diffx + diffy) - (advx + advy); + + if (stage == 1u) { + accum[scalar_idx] = rhs; + stage_out[scalar_idx] = base[scalar_idx] + p.alpha * rhs; + } else if (stage == 2u) { + accum[scalar_idx] = accum[scalar_idx] + p.weight * rhs; + stage_out[scalar_idx] = base[scalar_idx] + p.alpha * rhs; + } else if (stage == 3u) { + accum[scalar_idx] = accum[scalar_idx] + p.weight * rhs; + stage_out[scalar_idx] = base[scalar_idx] + p.alpha * rhs; + } else { + accum[scalar_idx] = accum[scalar_idx] + rhs; + base[scalar_idx] = base[scalar_idx] + p.dt_over_6 * accum[scalar_idx]; + stage_out[scalar_idx] = base[scalar_idx]; + } +} + +kernel void compute_fluxes_2d(device const float *stage [[buffer(0)]], + device const float *vx [[buffer(1)]], + device const float *vy [[buffer(2)]], + device float *diff_x [[buffer(3)]], + device float *adv_x [[buffer(4)]], + device float *diff_y [[buffer(5)]], + device float *adv_y [[buffer(6)]], + constant fa_params_t &p [[buffer(7)]], + uint2 gid [[thread_position_in_grid]]) +{ + fa_compute_fluxes_body(stage, vx, vy, diff_x, adv_x, diff_y, adv_y, p, + gid.x, gid.y); +} + +kernel void compute_fluxes_linear(device const float *stage [[buffer(0)]], + device const float *vx [[buffer(1)]], + device const float *vy [[buffer(2)]], + device float *diff_x [[buffer(3)]], + device float *adv_x [[buffer(4)]], + device float *diff_y [[buffer(5)]], + device float *adv_y [[buffer(6)]], + constant fa_params_t &p [[buffer(7)]], + uint gid [[thread_position_in_grid]]) +{ + const int flat = int(gid); + const int sx = int(p.nx + 2u); + const uint i = uint(flat - (flat / sx) * sx); + const uint j = uint(flat / sx); + fa_compute_fluxes_body(stage, vx, vy, diff_x, adv_x, diff_y, adv_y, p, i, j); +} + +kernel void update_stage_2d(device float *base [[buffer(0)]], + device const float *stage_in [[buffer(1)]], + device float *stage_out [[buffer(2)]], + device float *accum [[buffer(3)]], + device const float *diff_x [[buffer(4)]], + device const float *adv_x [[buffer(5)]], + device const float *diff_y [[buffer(6)]], + device const float *adv_y [[buffer(7)]], + constant fa_params_t &p [[buffer(8)]], + uint2 gid [[thread_position_in_grid]]) +{ + fa_update_stage_body(base, stage_in, stage_out, accum, diff_x, adv_x, diff_y, + adv_y, p, gid.x, gid.y); +} + +kernel void update_stage_linear(device float *base [[buffer(0)]], + device const float *stage_in [[buffer(1)]], + device float *stage_out [[buffer(2)]], + device float *accum [[buffer(3)]], + device const float *diff_x [[buffer(4)]], + device const float *adv_x [[buffer(5)]], + device const float *diff_y [[buffer(6)]], + device const float *adv_y [[buffer(7)]], + constant fa_params_t &p [[buffer(8)]], + uint gid [[thread_position_in_grid]]) +{ + const int flat = int(gid); + const int sx = int(p.nx + 2u); + const uint i = uint(flat - (flat / sx) * sx); + const uint j = uint(flat / sx); + fa_update_stage_body(base, stage_in, stage_out, accum, diff_x, adv_x, diff_y, + adv_y, p, i, j); +} diff --git a/benchmarks/2d-advection-metal/run.sh b/benchmarks/2d-advection-metal/run.sh new file mode 100755 index 0000000..fdeb0f4 --- /dev/null +++ b/benchmarks/2d-advection-metal/run.sh @@ -0,0 +1,69 @@ +#!/usr/bin/env bash +set -euo pipefail + +here="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +formal_root="$(cd "$here/../.." && pwd)" +lfortran_root="${LFORTRAN_ROOT:-"$formal_root/../lfortran"}" + +if [[ ! -f "$lfortran_root/pixi.toml" ]]; then + echo "LFortran checkout not found at $lfortran_root; set LFORTRAN_ROOT to its absolute path." >&2 + exit 2 +fi +lfortran_root="$(cd "$lfortran_root" && pwd)" +if [[ "${PIXI_PROJECT_ROOT:-}" != "$lfortran_root" || + "${PIXI_ENVIRONMENT_NAME:-}" != llvm11 ]]; then + exec pixi run --manifest-path "$lfortran_root/pixi.toml" -e llvm11 \ + bash "$here/run.sh" "$@" +fi + +if [[ "$(uname -s)" != Darwin || "$(uname -m)" != arm64 ]]; then + echo "the native Metal benchmark requires an Apple Silicon Mac" >&2 + exit 1 +fi + +gpu_guard="$lfortran_root/examples/nbody/gpu_run.sh" +if [[ ! -f "$gpu_guard" ]]; then + echo "LFortran GPU guard not found at $gpu_guard; use a complete LFortran checkout." >&2 + exit 2 +fi + +build_dir="$here/build" +binary="$build_dir/formal_advection" +mkdir -p "$build_dir" + +if [[ -n "${CC:-}" ]]; then + sdk_clang="$CC" +else + sdk_clang="$(xcrun --sdk macosx --find clang)" +fi +sdkroot="${SDKROOT:-$(xcrun --sdk macosx --show-sdk-path)}" +"$sdk_clang" -Wall -Wextra -Werror -O2 -fobjc-arc -isysroot "$sdkroot" \ + -framework Foundation -framework Metal \ + "$here/formal_advection.m" -o "$binary" + +if [[ $# -eq 0 ]]; then + set -- --mode cpu --cells 11,11 --steps 2 --dt 1e-4 +fi + +mode=cpu +for ((i = 1; i <= $#; ++i)); do + if [[ "${!i}" == "--mode" && $((i + 1)) -le $# ]]; then + j=$((i + 1)) + mode="${!j}" + break + fi +done +if [[ "$mode" != cpu && "$mode" != metal ]]; then + echo "--mode must be cpu or metal" >&2 + exit 2 +fi + +cd "$here" +if [[ "$mode" == metal ]]; then + if [[ "${FA_ALLOW_METAL:-0}" != 1 ]]; then + echo "Metal mode is guarded; set FA_ALLOW_METAL=1 to proceed." >&2 + exit 1 + fi + exec bash "$gpu_guard" "$binary" "$@" +fi +exec "$binary" "$@" diff --git a/benchmarks/2d-advection-metal/run_formal.sh b/benchmarks/2d-advection-metal/run_formal.sh new file mode 100755 index 0000000..1ebc8ed --- /dev/null +++ b/benchmarks/2d-advection-metal/run_formal.sh @@ -0,0 +1,50 @@ +#!/usr/bin/env bash +set -euo pipefail + +here="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +formal_root="$(cd "$here/../.." && pwd)" +lfortran_root="${LFORTRAN_ROOT:-"$formal_root/../lfortran"}" + +if [[ ! -f "$lfortran_root/pixi.toml" ]]; then + echo "LFortran checkout not found at $lfortran_root; set LFORTRAN_ROOT to its absolute path." >&2 + exit 2 +fi +lfortran_root="$(cd "$lfortran_root" && pwd)" +if [[ "${PIXI_PROJECT_ROOT:-}" != "$lfortran_root" || + "${PIXI_ENVIRONMENT_NAME:-}" != llvm11 ]]; then + exec pixi run --manifest-path "$lfortran_root/pixi.toml" -e llvm11 \ + bash "$here/run_formal.sh" "$@" +fi + +compiler="${LFORTRAN_COMPILER:-${LFORTRAN_BUILD_DIR:-"$lfortran_root/build/llvm11"}/src/bin/lfortran}" +gpu_guard="$lfortran_root/examples/nbody/gpu_run.sh" +if [[ ! -x "$compiler" ]]; then + echo "LFortran compiler not found at $compiler; build it with the llvm11 Pixi environment." >&2 + exit 2 +fi +if [[ ! -f "$formal_root/fpm.toml" || ! -f "$gpu_guard" ]]; then + echo "Formal checkout or LFortran GPU guard is incomplete; set LFORTRAN_ROOT explicitly." >&2 + exit 2 +fi + +compiler_hash="$(shasum -a 256 "$compiler" | awk '{print $1}')" +build_dir="$here/build/formal/compiler-${compiler_hash:0:16}" +mkdir -p "$build_dir" + +cd "$formal_root" +fpm run --compiler "$compiler" --build-dir "$build_dir" \ + --flag "--cpp --realloc-lhs-arrays --separate-compilation --gpu=metal" \ + --c-flag "-DFORMAL_ENABLE_GPU_BENCHMARK -I$lfortran_root/src/libasr/runtime" \ + --example 2D-advection-diffusion --runner cp -- "$build_dir/formal.metal" +if ! otool -L "$build_dir/formal.metal" | grep -q '/Metal.framework/'; then + echo "FPM did not produce a Metal-linked executable." >&2 + exit 1 +fi + +if [[ $# -eq 0 ]]; then + set -- --cells 2048,2048 --steps 1000 --no-csv +fi +export LFORTRAN_GPU_PROFILE=1 +GPU_TIMEOUT="${GPU_TIMEOUT:-30}" \ + bash "$gpu_guard" "$build_dir/formal.metal" --gpu-profile "$@" +printf 'Artifacts: %s\n' "$build_dir" From d3717125bb12e9b8641ea7a45cbaf49ccb9f4ac5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 20:30:47 -0600 Subject: [PATCH 06/18] Harden native Metal mode and working-set guards (cherry picked from commit 481b33e04feb1f9cd71c0c888b0fbbdf0e191f07) --- .../2d-advection-metal/formal_advection.m | 42 ++++++++++++++++--- benchmarks/2d-advection-metal/run.sh | 13 +++++- 2 files changed, 48 insertions(+), 7 deletions(-) diff --git a/benchmarks/2d-advection-metal/formal_advection.m b/benchmarks/2d-advection-metal/formal_advection.m index a2d79a6..13083cc 100644 --- a/benchmarks/2d-advection-metal/formal_advection.m +++ b/benchmarks/2d-advection-metal/formal_advection.m @@ -115,6 +115,14 @@ static size_t fa_checked_mul_size(size_t a, size_t b, const char *what) return a * b; } +static size_t fa_checked_add_size(size_t a, size_t b, const char *what) +{ + if (b > SIZE_MAX - a) { + fa_dief("%s size overflow", what); + } + return a + b; +} + static const char *fa_center_coord_text(size_t i, uint32_t n) { if (n == 51 && i < sizeof(FA_FORMAL_51_COORDS) / sizeof(FA_FORMAL_51_COORDS[0])) { @@ -288,11 +296,35 @@ static void fa_validate_working_set(const fa_config_t *cfg) (size_t)cfg->ny + 2, "x-face field"); const size_t yface_len = fa_checked_mul_size((size_t)cfg->nx + 2, (size_t)cfg->ny + 1, "y-face field"); - const size_t total_bytes = 4 * scalar_len * sizeof(float) + - 2 * xface_len * sizeof(float) + - 2 * yface_len * sizeof(float); - if (total_bytes > 512u * 1024u * 1024u) { - fa_die("estimated GPU working set exceeds 512 MiB"); + + const size_t scalar_bytes = fa_checked_mul_size(scalar_len, sizeof(float), + "scalar field bytes"); + const size_t xface_bytes = fa_checked_mul_size(xface_len, sizeof(float), + "x-face field bytes"); + const size_t yface_bytes = fa_checked_mul_size(yface_len, sizeof(float), + "y-face field bytes"); + const size_t host_fields_bytes = fa_checked_add_size( + fa_checked_add_size( + fa_checked_mul_size(scalar_bytes, 4, "host scalar buffers"), + xface_bytes, "host fields"), + yface_bytes, "host fields"); + const size_t warmup_backups_bytes = fa_checked_mul_size(scalar_bytes, 4, + "warmup backups"); + const size_t metal_buffers_bytes = fa_checked_add_size( + fa_checked_add_size( + fa_checked_mul_size(scalar_bytes, 4, "metal scalar buffers"), + fa_checked_mul_size(xface_bytes, 3, "metal x-face buffers"), + "metal buffers"), + fa_checked_mul_size(yface_bytes, 3, "metal y-face buffers"), + "metal buffers"); + const size_t total_bytes = fa_checked_add_size( + host_fields_bytes, + fa_checked_add_size(warmup_backups_bytes, metal_buffers_bytes, + "working set"), + "working set"); + const size_t max_bytes = (size_t)512u * 1024u * 1024u; + if (total_bytes > max_bytes) { + fa_die("estimated working set exceeds 512 MiB"); } } diff --git a/benchmarks/2d-advection-metal/run.sh b/benchmarks/2d-advection-metal/run.sh index fdeb0f4..b170abf 100755 --- a/benchmarks/2d-advection-metal/run.sh +++ b/benchmarks/2d-advection-metal/run.sh @@ -46,11 +46,20 @@ if [[ $# -eq 0 ]]; then fi mode=cpu +mode_seen=0 for ((i = 1; i <= $#; ++i)); do - if [[ "${!i}" == "--mode" && $((i + 1)) -le $# ]]; then + if [[ "${!i}" == "--mode" ]]; then + if ((i == $#)); then + echo "--mode must be followed by cpu or metal" >&2 + exit 2 + fi + if ((mode_seen)); then + echo "--mode may be specified only once" >&2 + exit 2 + fi j=$((i + 1)) mode="${!j}" - break + mode_seen=1 fi done if [[ "$mode" != cpu && "$mode" != metal ]]; then From 45435de07c9ac030069629031416e81732db47e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 20:38:03 -0600 Subject: [PATCH 07/18] docs: correct native advection provenance (cherry picked from commit 96c841a3c2c2aabbaec1588d0ed776a1ade0a1fc) --- benchmarks/2d-advection-metal/README.md | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/benchmarks/2d-advection-metal/README.md b/benchmarks/2d-advection-metal/README.md index 4a9832e..919f9fa 100644 --- a/benchmarks/2d-advection-metal/README.md +++ b/benchmarks/2d-advection-metal/README.md @@ -68,14 +68,27 @@ uses the same guarded launch path. It does not add native sources to FPM. ## Provenance and licensing -`formal_advection.m` and `formal_advection.metal` were copied verbatim from -the LFortran repository at commit +The numerical algorithm, FP32 CPU oracle, Metal kernels, two-deep scheduling, +and signed-index experiment come from the LFortran repository at commit `3d4c7941c5a153adb996ed62309225e871c06d9b` (`examples/formal_advection/formal_advection.m` and -`formal_advection.metal`). That commit is the source of the verified CPU -oracle, Metal kernels, two-deep scheduling, and signed-index experiment. -The adapted wrappers and this documentation are the Formal publication -integration. The upstream LFortran BSD terms are retained in +`examples/formal_advection/formal_advection.metal`). The provenance of each +native source file is recorded separately: + +- `formal_advection.metal` is byte-for-byte identical to the LFortran shader. + Its original and current Git blob hashes are both + `0e492038ab56e7fbb70789b724e934093392c3b0`. +- `formal_advection.m` retains that source's numerical method and CPU oracle, + but is not byte-for-byte identical. The Formal copy adds checked + `size_t` addition/multiplication and a 512 MiB native working-set guard + covering the host fields, warm-up backups, and Metal buffers. Its original + LFortran blob hash is `a2d79a6af54df22e4af9c4bffa85d5712942a2c8`; the + current Formal destination blob hash is + `13083cc9ab9547064b0d668552d747c1b8b203c2`. + +These C-side safety additions do not alter the numerical method or benchmark +claims. The adapted wrappers and this documentation are the Formal +publication integration. The upstream LFortran BSD terms are retained in [`LFortran-NOTICE.txt`](./LFortran-NOTICE.txt); the surrounding Formal tree remains covered by [`LICENSE.txt`](../../LICENSE.txt). From 1f7596e8be39b8dbcf981f2c5988cb54e97547a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 21:20:15 -0600 Subject: [PATCH 08/18] Require Metal opt-in and document paired native benchmarks (cherry picked from commit 846ffc5f0dc123411f9a685c569d2759f9a8a060) --- README.md | 10 +- benchmarks/2d-advection-metal/PERFORMANCE.md | 97 ++++++++++++++++++-- benchmarks/2d-advection-metal/README.md | 11 ++- benchmarks/2d-advection-metal/run_formal.sh | 9 ++ 4 files changed, 109 insertions(+), 18 deletions(-) diff --git a/README.md b/README.md index d74a94d..ae8fb1a 100644 --- a/README.md +++ b/README.md @@ -92,12 +92,12 @@ Formal GPU example through its wrappers: export LFORTRAN_ROOT=/path/to/lfortran bash benchmarks/2d-advection-metal/run.sh # or: -bash benchmarks/2d-advection-metal/run_formal.sh +FA_ALLOW_METAL=1 bash benchmarks/2d-advection-metal/run_formal.sh ``` -The native wrapper defaults to a CPU smoke case; Metal is explicitly guarded -and the measured 2048 x 2048, 1000-step case is selected by passing its -arguments. See the benchmark README for the exact provenance, source guard, -and reproducibility contract. +The native wrapper defaults to a CPU smoke case; the Formal wrapper defaults +to the 2048 x 2048, 1000-step GPU case. Both Metal entry points require an +Apple Silicon Mac and explicit `FA_ALLOW_METAL=1` permission. See the benchmark +README for the exact provenance, source guard, and reproducibility contract. ### Other example programs See this repository's [example](./example) subdirectory for additional demonstrations of diff --git a/benchmarks/2d-advection-metal/PERFORMANCE.md b/benchmarks/2d-advection-metal/PERFORMANCE.md index 5aad44f..0b33843 100644 --- a/benchmarks/2d-advection-metal/PERFORMANCE.md +++ b/benchmarks/2d-advection-metal/PERFORMANCE.md @@ -24,22 +24,99 @@ implementations. On the Apple M4 Max used for the original comparison, the measured Gaussian case was 2048 x 2048 cells for 1000 steps with the example's scaled timestep -(approximately `6.20126741e-8`). Run the native reference with: +(approximately `6.20126741e-8`). The native variant used direct 2D indexing, +**128 x 1 threadgroups**, and the **same raw Formal input**, not its default +32 x 4 threadgroups or independently generated analytic input. + +Run the following from the Formal repository root on an Apple Silicon Mac, +in one Bash session. Set `LFORTRAN_ROOT` to an absolute checkout with its +LLVM 11 Pixi compiler already built. Metal permission is deliberately +explicit; neither wrapper enables it automatically. ```bash -LFORTRAN_ROOT=/path/to/lfortran \ +set -euo pipefail +export LFORTRAN_ROOT=/path/to/lfortran +export FA_ALLOW_METAL=1 GPU_TIMEOUT=30 GPU_MAXMEM_MB=4096 +nx=2048 +steps=1000 +pair_dir="$PWD/benchmarks/2d-advection-metal/build/paired-${nx}-${steps}" +mkdir -p "$pair_dir" + +pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ + bash benchmarks/2d-advection-metal/run_formal.sh \ + --solver workspace --order 4 --cells "$nx,$nx" --steps "$steps" --no-csv \ + --dump-input "$pair_dir/input.fadv2d" \ + --dump-final-raw "$pair_dir/formal-final.fadv2d" \ + > "$pair_dir/formal.log" 2>&1 + +pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ bash benchmarks/2d-advection-metal/run.sh \ - --mode metal --cells 2048,2048 --steps 1000 + --mode metal --input "$pair_dir/input.fadv2d" \ + --index-mode 2d --tg 128,1 --dump-bin "$pair_dir/native-final.bin" \ + > "$pair_dir/native.log" 2>&1 + +cat "$pair_dir/formal.log" "$pair_dir/native.log" +``` + +Formal exports the initial fields before warmup or solve timing. Native +`--input` reads those fields and their exact grid, step count, and FP32 +timestep from the snapshot; do not regenerate the input or round the +timestep. Absolute artifact paths are important: the Formal wrapper runs +from the repository root, while the native wrapper runs from this directory. +Both wrappers pass the actual executable directly to LFortran's existing +GPU guard. The guard is required for every Metal run and serializes users +through `/tmp/gpurun.lock.d`. +Check that Formal reports `gpu_timing_available=T` and positive +`gpu_busy_seconds`, and native reports `mode=metal`, `index_mode=2d`, +`tg=128,1`, and positive `gpu_busy`. These solve timings exclude compilation, +initialization, warmup/JIT, and field I/O. Neither command writes CSV. + +For a bounded command-wiring check, change only `nx=17` and `steps=9` above. +That is a reduced check, not a reproduction of the historical timings. + +To compare the complete FP32 final fields, including fixed boundaries, use +the same shell variables. Formal's final file has a three-line FADV2D1 +header; native's binary has only the scalar payload. This check rejects +incorrect element counts and non-finite fields before computing errors: + +```bash +pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ + python - "$pair_dir" "$nx" "$steps" > "$pair_dir/compare.log" 2>&1 <<'PY' +import math +from pathlib import Path +import sys +import numpy as np + +pair_dir = Path(sys.argv[1]) +nx, steps = map(int, sys.argv[2:]) +with (pair_dir / "formal-final.fadv2d").open("rb") as stream: + assert stream.readline().strip() == b"FADV2D1", "invalid Formal header" + x, y, n, dt = stream.readline().split() + assert (int(x), int(y), int(n)) == (nx, nx, steps), "configuration mismatch" + assert math.isfinite(float(dt)) and float(dt) > 0, "invalid timestep" + assert stream.readline().strip() == b"mode raw-final", "invalid final mode" + formal = np.frombuffer(stream.read(), dtype="`, `--dump-bin [path]`, flattened-index experiment. The verified comparison found it slower on the small cache-resident case; no claim that unsigned indexing is faster is made. -To run the actual Formal GPU example rather than the native oracle: +To run the actual Formal GPU example rather than the native oracle, the same +Apple Silicon and explicit Metal opt-in requirements apply: ```bash -LFORTRAN_ROOT=/path/to/lfortran \ +FA_ALLOW_METAL=1 LFORTRAN_ROOT=/path/to/lfortran \ bash benchmarks/2d-advection-metal/run_formal.sh ``` -That wrapper compiles the existing +That wrapper rejects missing permission or an unsupported platform before +compilation. It defaults to 2048 x 2048 cells and 1000 steps, compiles `example/2D-advection-diffusion.F90` with LFortran, checks Metal linkage, and uses the same guarded launch path. It does not add native sources to FPM. +See [`PERFORMANCE.md`](./PERFORMANCE.md#reproduction) for the paired +shared-input recipe and the explicit native launch geometry used in the +historical comparison. ## Provenance and licensing diff --git a/benchmarks/2d-advection-metal/run_formal.sh b/benchmarks/2d-advection-metal/run_formal.sh index 1ebc8ed..b11b375 100755 --- a/benchmarks/2d-advection-metal/run_formal.sh +++ b/benchmarks/2d-advection-metal/run_formal.sh @@ -16,6 +16,15 @@ if [[ "${PIXI_PROJECT_ROOT:-}" != "$lfortran_root" || bash "$here/run_formal.sh" "$@" fi +if [[ "$(uname -s)" != Darwin || "$(uname -m)" != arm64 ]]; then + echo "the Formal Metal benchmark requires an Apple Silicon Mac" >&2 + exit 1 +fi +if [[ "${FA_ALLOW_METAL:-0}" != 1 ]]; then + echo "Metal mode is guarded; set FA_ALLOW_METAL=1 to proceed." >&2 + exit 1 +fi + compiler="${LFORTRAN_COMPILER:-${LFORTRAN_BUILD_DIR:-"$lfortran_root/build/llvm11"}/src/bin/lfortran}" gpu_guard="$lfortran_root/examples/nbody/gpu_run.sh" if [[ ! -x "$compiler" ]]; then From 4d6a2a1716f73c69e279c77c570de3057c18fa1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 21:56:26 -0600 Subject: [PATCH 09/18] Reject unevolved native benchmark field comparisons (cherry picked from commit 5a3bc758d8a915d3864e6ae667a4f3625dea27bb) --- benchmarks/2d-advection-metal/PERFORMANCE.md | 75 +++--- example/scripts/compare_2d_fields.py | 9 +- example/scripts/compare_2d_raw_fields.py | 230 ++++++++++++++++++ .../validate_2d_advection_workspace.py | 26 +- 4 files changed, 301 insertions(+), 39 deletions(-) create mode 100644 example/scripts/compare_2d_raw_fields.py diff --git a/benchmarks/2d-advection-metal/PERFORMANCE.md b/benchmarks/2d-advection-metal/PERFORMANCE.md index 0b33843..86f959f 100644 --- a/benchmarks/2d-advection-metal/PERFORMANCE.md +++ b/benchmarks/2d-advection-metal/PERFORMANCE.md @@ -38,13 +38,14 @@ set -euo pipefail export LFORTRAN_ROOT=/path/to/lfortran export FA_ALLOW_METAL=1 GPU_TIMEOUT=30 GPU_MAXMEM_MB=4096 nx=2048 +ny=2048 steps=1000 -pair_dir="$PWD/benchmarks/2d-advection-metal/build/paired-${nx}-${steps}" +pair_dir="$PWD/benchmarks/2d-advection-metal/build/paired-${nx}x${ny}-${steps}" mkdir -p "$pair_dir" pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ bash benchmarks/2d-advection-metal/run_formal.sh \ - --solver workspace --order 4 --cells "$nx,$nx" --steps "$steps" --no-csv \ + --solver workspace --order 4 --cells "$nx,$ny" --steps "$steps" --no-csv \ --dump-input "$pair_dir/input.fadv2d" \ --dump-final-raw "$pair_dir/formal-final.fadv2d" \ > "$pair_dir/formal.log" 2>&1 @@ -71,47 +72,51 @@ Check that Formal reports `gpu_timing_available=T` and positive `tg=128,1`, and positive `gpu_busy`. These solve timings exclude compilation, initialization, warmup/JIT, and field I/O. Neither command writes CSV. -For a bounded command-wiring check, change only `nx=17` and `steps=9` above. +For a bounded non-square command-wiring check, use `nx=17`, `ny=13`, and +`steps=9` above. That is a reduced check, not a reproduction of the historical timings. -To compare the complete FP32 final fields, including fixed boundaries, use -the same shell variables. Formal's final file has a three-line FADV2D1 -header; native's binary has only the scalar payload. This check rejects -incorrect element counts and non-finite fields before computing errors: +To validate the complete FP32 fields, use the same shell variables and the +shared raw/field comparator below. It reads the initial scalar and velocity +payloads as well as both finals, checks their sizes and finiteness, and +requires matching grid, step count, and FP32 timestep metadata. Fixed scalar +boundaries must equal the initial values exactly. Final-field agreement uses +the workspace validator's tolerances (`rtol=1e-6`, `atol=1e-8`) at every point, +including boundaries, and reports the relative L2 error. + +For this nonstationary Gaussian, **keep `--evolution`**: both final fields +must be nonzero and must differ meaningfully from the shared initial field. +The maximum change must exceed `atol + rtol * max(abs(initial))`, and the +relative L2 change must exceed `rtol`. This rejects unchanged initial fields, +zeroed finals, and roundoff-only changes even when the two finals agree. +The option is deliberately not a requirement for legitimate stationary or +zero-state tests. Positive GPU time establishes execution, not numerical +correctness; both telemetry and field validation are required. ```bash pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ - python - "$pair_dir" "$nx" "$steps" > "$pair_dir/compare.log" 2>&1 <<'PY' -import math -from pathlib import Path -import sys -import numpy as np - -pair_dir = Path(sys.argv[1]) -nx, steps = map(int, sys.argv[2:]) -with (pair_dir / "formal-final.fadv2d").open("rb") as stream: - assert stream.readline().strip() == b"FADV2D1", "invalid Formal header" - x, y, n, dt = stream.readline().split() - assert (int(x), int(y), int(n)) == (nx, nx, steps), "configuration mismatch" - assert math.isfinite(float(dt)) and float(dt) > 0, "invalid timestep" - assert stream.readline().strip() == b"mode raw-final", "invalid final mode" - formal = np.frombuffer(stream.read(), dtype=" "$pair_dir/compare.log" 2>&1 cat "$pair_dir/compare.log" ``` +For small cases, `--cpu-final ` additionally checks both +finals against the retained native CPU oracle, with the same tolerances, +evolution requirement, and fixed-boundary checks. Comparator regressions can +be checked without a GPU: + +```bash +mkdir -p benchmarks/2d-advection-metal/build +pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ + python example/scripts/compare_2d_raw_fields.py --self-test \ + > benchmarks/2d-advection-metal/build/compare-selftest.log 2>&1 +cat benchmarks/2d-advection-metal/build/compare-selftest.log +``` + +These synthetic comparator tests are not solver validation. + ## Historical paired results The original measurement used three interleaved pairs with the configuration diff --git a/example/scripts/compare_2d_fields.py b/example/scripts/compare_2d_fields.py index 39e6253..ecc8ab1 100755 --- a/example/scripts/compare_2d_fields.py +++ b/example/scripts/compare_2d_fields.py @@ -122,8 +122,7 @@ def load_and_validate(reference_path: Path, observed_path: Path): return reference, observed -def run_comparison(reference_path: Path, observed_path: Path, rtol: float, atol: float, evolution: bool): - reference, observed = load_and_validate(reference_path, observed_path) +def report_comparison(reference, observed, rtol: float, atol: float, evolution: bool): if evolution: max_abs_diff, max_rel_diff, relative_l2, changed = summarize_field_differences(reference, observed) else: @@ -134,6 +133,12 @@ def run_comparison(reference_path: Path, observed_path: Path, rtol: float, atol: print(f"changed={str(changed).lower()}") if evolution and not changed: raise SystemExit("field did not change") + return max_abs_diff, max_rel_diff, relative_l2, changed + + +def run_comparison(reference_path: Path, observed_path: Path, rtol: float, atol: float, evolution: bool): + reference, observed = load_and_validate(reference_path, observed_path) + report_comparison(reference, observed, rtol, atol, evolution) def write_csv(path: Path, xs, ys, func): diff --git a/example/scripts/compare_2d_raw_fields.py b/example/scripts/compare_2d_raw_fields.py new file mode 100644 index 0000000..abf7c9c --- /dev/null +++ b/example/scripts/compare_2d_raw_fields.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import math +import sys +from pathlib import Path + +from compare_2d_fields import report_comparison +from validate_2d_advection_workspace import parse_raw_headered, parse_raw_headerless + + +def check_boundaries(initial, final, label): + nx, ny = len(initial) - 2, len(initial[0]) - 2 + boundary = [(i, j) for i in range(nx + 2) for j in (0, ny + 1)] + boundary += [(i, j) for i in (0, nx + 1) for j in range(1, ny + 1)] + for i, j in boundary: + if final[i][j] != initial[i][j]: + raise SystemExit(f"{label}: fixed boundary changed at ({i},{j})") + + +def compare_raw_fields(initial_path, formal_path, native_path, cells, steps, + rtol=1e-6, atol=1e-8, evolution=False, cpu_path=None): + if len(cells) != 2 or min(cells) < 9 or steps <= 0: + raise SystemExit("expected order-4 cells >= 9 and a positive step count") + if not all(math.isfinite(t) and t >= 0 for t in (rtol, atol)): + raise SystemExit("--rtol and --atol must be finite and non-negative") + initial_snapshot = parse_raw_headered(initial_path, expected_mode="raw") + final_snapshot = parse_raw_headered(formal_path, expected_mode="raw-final") + nx, ny, n, dt, initial = initial_snapshot + if (nx, ny, n) != (*cells, steps): + raise SystemExit("initial configuration mismatch") + if final_snapshot[:4] != initial_snapshot[:4]: + raise SystemExit("Formal final configuration differs from shared input") + formal = final_snapshot[4] + native = parse_raw_headerless(native_path, nx, ny) + finals = [("Formal", formal), ("native", native)] + if cpu_path is not None: + finals.append(("CPU oracle", parse_raw_headerless(cpu_path, nx, ny))) + + initial_scale = max(abs(value) for row in initial for value in row) + if evolution and initial_scale == 0: + raise SystemExit("--evolution requires a nonzero initial field") + print(f"cells={nx},{ny} steps={n} dt={dt:.16e} elements={(nx + 2) * (ny + 2)}") + for label, final in finals: + if evolution and not any(value != 0 for row in final for value in row): + raise SystemExit(f"{label}: zero final field from nonzero input") + check_boundaries(initial, final, label) + if evolution: + print(f"{label} evolution:") + max_abs, _, relative_l2, _ = report_comparison(initial, final, rtol, atol, True) + if max_abs <= atol + rtol * initial_scale or relative_l2 <= rtol: + raise SystemExit(f"{label}: evolution does not exceed comparison tolerances") + print("Formal versus native:") + report_comparison(formal, native, rtol, atol, False) + if cpu_path is not None: + for label, final in finals[:2]: + print(f"CPU oracle versus {label}:") + report_comparison(finals[2][1], final, rtol, atol, False) + print("full-field comparison passed; fixed boundaries preserved") + + +def self_test(): + import contextlib + import io + import shutil + import struct + + scratch = Path(__file__).resolve().parents[2] / "build/formal-gpu/compare_2d_raw_fields-selftest" + scratch.mkdir(parents=True, exist_ok=False) + nx, ny, steps = 13, 11, 3 + count = (nx + 2) * (ny + 2) + scalar = [0.5 + math.exp(-((i - 5) ** 2 + (j - 4) ** 2) / 8) + for j in range(ny + 2) for i in range(nx + 2)] + evolved = scalar.copy() + for j in range(1, ny + 1): + for i in range(1, nx + 1): + evolved[i + j * (nx + 2)] *= 0.99 + velocities = [0.1] * ((nx + 1) * (ny + 2) + (nx + 2) * (ny + 1)) + input_header = f"FADV2D1\n{nx} {ny} {steps} 0.0001\nmode raw\n".encode() + final_header = input_header.replace(b"mode raw\n", b"mode raw-final\n") + + def pack(values): + return struct.pack(f"<{len(values)}f", *values) + + initial_bytes = input_header + pack(scalar + velocities) + formal_bytes = final_header + pack(evolved) + native_bytes = pack(evolved) + paths = [scratch / name for name in ("input.fadv2d", "formal.fadv2d", "native.bin")] + checks = 0 + + def check(label, replacements=None, error=None, **kwargs): + nonlocal checks + data = [initial_bytes, formal_bytes, native_bytes] + for index, payload in (replacements or {}).items(): + data[index] = payload + for path, payload in zip(paths, data): + if payload is None: + path.unlink(missing_ok=True) + else: + path.write_bytes(payload) + options = dict(cells=(nx, ny), steps=steps, evolution=True) + options.update(kwargs) + try: + with contextlib.redirect_stdout(io.StringIO()): + compare_raw_fields(*paths, **options) + except (SystemExit, OSError, ValueError) as exc: + if error is None or error not in str(exc): + raise SystemExit(f"self-test {label}: unexpected failure: {exc}") from exc + else: + if error is not None: + raise SystemExit(f"self-test {label}: expected rejection") + checks += 1 + print(f"PASS {label}") + + try: + check("synthetic evolved pair") + check("unchanged initial", {1: final_header + pack(scalar), 2: pack(scalar)}, + "field did not change") + zeros = pack([0.0] * count) + check("both finals zero", {1: final_header + zeros, 2: zeros}, "zero final field") + check("Formal unchanged", {1: final_header + pack(scalar)}, "field did not change") + check("native unchanged", {2: pack(scalar)}, "field did not change") + tiny = scalar.copy() + tiny[nx + 3] += 1e-7 + check("roundoff is not evolution", {1: final_header + pack(tiny), 2: pack(tiny)}, + "evolution does not exceed") + zero_input = input_header + zeros + pack(velocities) + check("zero initial cannot require evolution", {0: zero_input}, "nonzero initial") + check("stationary zero remains valid without evolution", + {0: zero_input, 1: final_header + zeros, 2: zeros}, evolution=False) + drift = evolved.copy() + drift[0] += 0.25 + check("both boundaries drift", {1: final_header + pack(drift), 2: pack(drift)}, + "fixed boundary changed") + check("native boundary drift", {2: pack(drift)}, "native: fixed boundary changed") + wrong = evolved.copy() + wrong[nx + 3] += 0.05 + check("finite interior mismatch", {2: pack(wrong)}, "field mismatch") + for index, name in enumerate(("initial", "Formal", "native")): + check(f"missing {name}", {index: None}, "No such file") + check(f"truncated {name}", {index: (initial_bytes, formal_bytes, native_bytes)[index][:-4]}, + "payload size") + check(f"extra {name}", {index: (initial_bytes, formal_bytes, native_bytes)[index] + b"\0" * 4}, + "payload size") + check("partial FP32 element", {2: native_bytes[:-1]}, "payload size") + for value in (math.nan, math.inf): + for index, name in enumerate(("initial scalar", "Formal scalar", "native scalar", "x velocity", "y velocity")): + if index in (0, 3, 4): + values = scalar + velocities + offset = 0 if index == 0 else count + if index == 4: + offset += (nx + 1) * (ny + 2) + values[offset] = value + change = {0: input_header + pack(values)} + else: + values = evolved.copy() + values[nx + 3] = value + change = {index: (final_header if index == 1 else b"") + pack(values)} + check(f"{name} {value}", change, "non-finite") + check("expected cells", error="initial configuration mismatch", cells=(ny, nx)) + check("expected steps", error="initial configuration mismatch", steps=4) + check("zero expected steps", error="positive step count", steps=0) + check("wrong Formal steps", {1: formal_bytes.replace(b"13 11 3", b"13 11 4")}, + "configuration differs") + check("wrong Formal timestep", {1: formal_bytes.replace(b"0.0001", b"0.0002")}, + "configuration differs") + for dt in (b"nan", b"inf", b"0", b"-1", b"1e-100"): + check(f"invalid initial dt {dt!r}", {0: initial_bytes.replace(b"0.0001", dt)}, + "timestep must be finite and positive") + check("overflowing initial timestep", {0: initial_bytes.replace(b"0.0001", b"1e100")}, + "invalid FP32 timestep") + rounded_dt = f"{struct.unpack(' Date: Sat, 17 Oct 2026 22:52:41 -0600 Subject: [PATCH 10/18] fix: pass a procedure pointer to workspace reconstruction Match the scalar_2D_t constructor's callback interface using the existing pure module initializer. Preserve the reconstruction buffer, its lifetime, and every field value. (cherry picked from commit a0e4e8d0387855c6d12df876e29c354eda7aebd2) --- example/formal_advection_workspace.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/example/formal_advection_workspace.F90 b/example/formal_advection_workspace.F90 index 69a0510..bbab840 100644 --- a/example/formal_advection_workspace.F90 +++ b/example/formal_advection_workspace.F90 @@ -1,6 +1,6 @@ module formal_advection_workspace use iso_fortran_env, only: real32, real64, int64 - use tensors_2D_m, only: scalar_2D_t + use tensors_2D_m, only: scalar_2D_t, scalar_2D_initializer_i implicit none private public :: fa_config_t, fa_options_t, fa_fields_t, fa_parse_args @@ -504,9 +504,11 @@ function fa_reconstruct_scalar(values, order, cells, x_min, x_max) result(scalar integer, intent(in) :: cells(:) real(real32), intent(in) :: x_min(:), x_max(:) type(scalar_2D_t) :: scalar + procedure(scalar_2D_initializer_i), pointer :: initializer + initializer => fa_scalar_reconstruction_initializer fa_scalar_reconstruction_values = values - scalar = scalar_2D_t(fa_scalar_reconstruction_initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) + scalar = scalar_2D_t(initializer, order=order, cells=cells, x_min=x_min, x_max=x_max) deallocate(fa_scalar_reconstruction_values) end function fa_reconstruct_scalar From 7393d21ef267fadf39145550a1e86b52854b2735 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sat, 17 Oct 2026 23:25:59 -0600 Subject: [PATCH 11/18] fix: export workspace array initialization and raw output (cherry picked from commit 72b4b2839272b1af32e9392453c96a4beaf59672) --- example/formal_advection_workspace.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/example/formal_advection_workspace.F90 b/example/formal_advection_workspace.F90 index bbab840..0bf5c0e 100644 --- a/example/formal_advection_workspace.F90 +++ b/example/formal_advection_workspace.F90 @@ -5,6 +5,7 @@ module formal_advection_workspace private public :: fa_config_t, fa_options_t, fa_fields_t, fa_parse_args public :: fa_init_fields, fa_advance_rk4, fa_write_csv, fa_reconstruct_scalar + public :: fa_init_from_arrays, fa_write_final_raw public :: fa_validate_cfg, fa_validate_working_set, fa_read_raw_input public :: fa_now_seconds From 7bb6c82229bc1932992ac3df0ece1752967882be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sun, 18 Oct 2026 00:35:33 -0600 Subject: [PATCH 12/18] docs: record legacy Fortran compiler limits (cherry picked from commit db4b32ae265280d74c97edecd7dda30b9796c1f8) --- benchmarks/2d-advection-metal/README.md | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/benchmarks/2d-advection-metal/README.md b/benchmarks/2d-advection-metal/README.md index 64d505c..4b3bf9c 100644 --- a/benchmarks/2d-advection-metal/README.md +++ b/benchmarks/2d-advection-metal/README.md @@ -71,6 +71,22 @@ See [`PERFORMANCE.md`](./PERFORMANCE.md#reproduction) for the paired shared-input recipe and the explicit native launch geometry used in the historical comparison. +## Known Fortran compiler limitations + +The standalone C/Metal reference is independent of the Fortran compiler. +The Fortran example has two known compiler-version failures in +[branch CI](https://github.com/certik/formal/actions/runs/39102014517): + +- CI's Flang 21 crashes while lowering `do concurrent ... local(rhs, total)`. +- LFortran 0.64 fails LLVM type verification for `numerical_arguments_t` + while compiling the example's argument and raw-input helpers. + +Both failures occur before the solver runs. Current LFortran builds and the +other completed Flang configurations compile the portable example; real GPU +comparisons verify evolved full fields and unchanged boundaries. These older +compiler limitations are recorded without changing the numerical method or +disabling their CI checks. + ## Provenance and licensing The numerical algorithm, FP32 CPU oracle, Metal kernels, two-deep scheduling, From 14da124fabe8c5bacf175fef5986043b2f934b4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sun, 18 Oct 2026 05:59:20 -0600 Subject: [PATCH 13/18] docs: clarify legacy macOS CI coverage (cherry picked from commit adf06916aa56805e53b35d39ab47c216cd0d02bb) --- benchmarks/2d-advection-metal/README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/benchmarks/2d-advection-metal/README.md b/benchmarks/2d-advection-metal/README.md index 4b3bf9c..59aa9c4 100644 --- a/benchmarks/2d-advection-metal/README.md +++ b/benchmarks/2d-advection-metal/README.md @@ -87,6 +87,11 @@ comparisons verify evolved full fields and unchanged boundaries. These older compiler limitations are recorded without changing the numerical method or disabling their CI checks. +That CI run also failed to build the Z3 dependency on macOS 14. The two Intel +macOS jobs were cancelled during LLVM installation, without a more specific +reported reason. None of these three jobs reached Formal compilation or tests; +they do not provide validation of those platforms. + ## Provenance and licensing The numerical algorithm, FP32 CPU oracle, Metal kernels, two-deep scheduling, From 70e73ee51f72ebaca16544d834bdfd98a9da66a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sun, 18 Oct 2026 09:15:15 -0600 Subject: [PATCH 14/18] Make advection benchmarks reproducible with managed FPM --- README.md | 4 + benchmarks/.gitignore | 1 + benchmarks/2d-advection-metal/PERFORMANCE.md | 181 ++++++++++++++++++- benchmarks/2d-advection-metal/README.md | 69 +++++-- benchmarks/2d-advection-metal/run.sh | 78 ++++---- benchmarks/2d-advection-metal/run_formal.sh | 35 ++-- benchmarks/pixi.lock | 99 ++++++++++ benchmarks/pixi.toml | 8 + benchmarks/toolchain.sh | 51 ++++++ 9 files changed, 436 insertions(+), 90 deletions(-) create mode 100644 benchmarks/.gitignore create mode 100644 benchmarks/pixi.lock create mode 100644 benchmarks/pixi.toml create mode 100644 benchmarks/toolchain.sh diff --git a/README.md b/README.md index ae8fb1a..10a3cf6 100644 --- a/README.md +++ b/README.md @@ -98,6 +98,10 @@ The native wrapper defaults to a CPU smoke case; the Formal wrapper defaults to the 2048 x 2048, 1000-step GPU case. Both Metal entry points require an Apple Silicon Mac and explicit `FA_ALLOW_METAL=1` permission. See the benchmark README for the exact provenance, source guard, and reproducibility contract. +The optional benchmark wrappers provision their own locked FPM 0.13.0 +environment from `benchmarks/pixi.toml`, reuse the selected LFortran +checkout's managed tools, and select the real macOS SDK compiler after Pixi +activation. Ordinary FPM users and the compiler support matrix are unchanged. ### Other example programs See this repository's [example](./example) subdirectory for additional demonstrations of diff --git a/benchmarks/.gitignore b/benchmarks/.gitignore new file mode 100644 index 0000000..604126f --- /dev/null +++ b/benchmarks/.gitignore @@ -0,0 +1 @@ +/.pixi/ diff --git a/benchmarks/2d-advection-metal/PERFORMANCE.md b/benchmarks/2d-advection-metal/PERFORMANCE.md index 86f959f..cec4a69 100644 --- a/benchmarks/2d-advection-metal/PERFORMANCE.md +++ b/benchmarks/2d-advection-metal/PERFORMANCE.md @@ -32,26 +32,29 @@ Run the following from the Formal repository root on an Apple Silicon Mac, in one Bash session. Set `LFORTRAN_ROOT` to an absolute checkout with its LLVM 11 Pixi compiler already built. Metal permission is deliberately explicit; neither wrapper enables it automatically. +The Formal wrapper provisions FPM 0.13.0 from the benchmark-local Pixi +manifest/lockfile; it does not assume FPM exists in LFortran's environment. +Both wrappers select the real SDK compiler **after** Pixi activation. ```bash set -euo pipefail export LFORTRAN_ROOT=/path/to/lfortran +export LFORTRAN_COMPILER="$LFORTRAN_ROOT/build/llvm11/src/bin/lfortran" export FA_ALLOW_METAL=1 GPU_TIMEOUT=30 GPU_MAXMEM_MB=4096 nx=2048 ny=2048 steps=1000 pair_dir="$PWD/benchmarks/2d-advection-metal/build/paired-${nx}x${ny}-${steps}" -mkdir -p "$pair_dir" +mkdir -p "$pair_dir/scratch" +export TMPDIR="$pair_dir/scratch" -pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ - bash benchmarks/2d-advection-metal/run_formal.sh \ +bash benchmarks/2d-advection-metal/run_formal.sh \ --solver workspace --order 4 --cells "$nx,$ny" --steps "$steps" --no-csv \ --dump-input "$pair_dir/input.fadv2d" \ --dump-final-raw "$pair_dir/formal-final.fadv2d" \ > "$pair_dir/formal.log" 2>&1 -pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ - bash benchmarks/2d-advection-metal/run.sh \ +bash benchmarks/2d-advection-metal/run.sh \ --mode metal --input "$pair_dir/input.fadv2d" \ --index-mode 2d --tg 128,1 --dump-bin "$pair_dir/native-final.bin" \ > "$pair_dir/native.log" 2>&1 @@ -95,7 +98,7 @@ correctness; both telemetry and field validation are required. ```bash pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ - python example/scripts/compare_2d_raw_fields.py \ + "$LFORTRAN_ROOT/.pixi/envs/llvm11/bin/python" example/scripts/compare_2d_raw_fields.py \ "$pair_dir/input.fadv2d" "$pair_dir/formal-final.fadv2d" \ "$pair_dir/native-final.bin" --cells "$nx,$ny" --steps "$steps" \ --rtol 1e-6 --atol 1e-8 --evolution > "$pair_dir/compare.log" 2>&1 @@ -110,13 +113,172 @@ be checked without a GPU: ```bash mkdir -p benchmarks/2d-advection-metal/build pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ - python example/scripts/compare_2d_raw_fields.py --self-test \ + "$LFORTRAN_ROOT/.pixi/envs/llvm11/bin/python" example/scripts/compare_2d_raw_fields.py --self-test \ > benchmarks/2d-advection-metal/build/compare-selftest.log 2>&1 cat benchmarks/2d-advection-metal/build/compare-selftest.log ``` These synthetic comparator tests are not solver validation. +## Controlled CPU measurements + +Measured on **18 October 2026, Apple M4 Max, 64 GiB**, macOS 15.8 +(Darwin 24.6), AC power, `powermode=0`. Timed processes used +`OMP_NUM_THREADS=OPENBLAS_NUM_THREADS=VECLIB_MAXIMUM_THREADS=1`; builds used +at most four jobs. This was a shared desktop, without core affinity or locked +clocks: load averages were about 2.4–2.5 during the sequential CPU run. +The CPU has 12 performance and four efficiency cores; no core class was forced. +Three samples are descriptive measurements, not a statistical significance +test or a claim of a universal maximum. + +The optimized numerical sources are unchanged from `adf06916`. +Two isolated **source-only** snapshots of the combined branch were used: +the control differs only in these files, restored to correctness prerequisite +`363f8999` (their contents also match the parent of `3866131`): + +```text +src/formal/differential_operators_1D_m.F90 +src/formal/divergence_operator_1D_s.F90 +src/formal/gradient_operator_1D_s.F90 +``` + +No compiled modules, build caches or environments were copied. Both used the +same instrumented example and original `--solver operators`; the workspace +comparison then used the **same folded executable** with `--solver workspace`. +All inputs within each compiler/workload were byte-identical, including +velocities, grid, step count and FP32 timestep. Raw `--input` is workspace-only, +so the operator controls used the same analytic parameters and compared their +exported raw inputs, rather than passing an unsupported option. + +Toolchains and matching flags: + +- **GFortran 16.2.0**, from the clean LFortran `llvm11` Pixi environment: + FPM `--profile release --flag "-fcoarray=single -B/usr/bin/ -ffree-line-length-none"`. + The release profile adds `-O3 -Wimplicit-interface -fPIC -fmax-errors=1 + -funroll-loops`; it does not enable `-ffast-math`. +- **Clean LFortran `2247d35082ea190408486b93c12ec2ead4771735`**, version + `0.67.0-922-g2247d35082`, LLVM 11.1.0: + `--flag "--cpp --realloc-lhs-arrays --separate-compilation --fast"`, + **without** FPM `--profile release`. FPM 0.13.0's LFortran release profile + emitted a literal `flag_lfortran_opt`, which the linker rejected. The + explicit `--fast` invocation avoids that tooling defect without changing + FPM or the compiler. Compiler SHA-256: + `8a1064336ee4922499cb81674c1fc04fbdd376a4b43a8388d993b8e775f22f8d`. + This is an unpatched Debug compiler built with `-g -O2` and assertions; + those compiler-build flags are distinct from the generated program's `--fast`. +- FPM **0.13.0** from Formal's benchmark lockfile; real Apple Clang 17 and + SDK 26.2 for C/linking; Python 3.12/NumPy 1.26.4 from clean LFortran's Pixi. + +One calibration at 51×51/500, `dt=1e-4`, selected the following **frozen** +step counts to keep solves around 0.25–2 s. Each invocation performs the +unchanged untimed one-step warmup and restores the initial state. Three +interleaved triplets used original/folded/workspace, then the reverse order, +then original/folded/workspace. Solve times exclude build, initial setup, +warmup and file I/O; the workspace's once-per-advance velocity scratch +allocation/copy is still inside the solve. + +| Compiler / workload | Original operators, median [min, max] s | Folded operators, median [min, max] s | Folded workspace, median [min, max] s | +|---|---:|---:|---:| +| GFortran, 51×51 / 2000 steps | 0.769787 [0.769689, 0.770542] | 0.769291 [0.762533, 0.787986] | 0.264657 [0.263707, 0.265147] | +| LFortran `--fast`, 51×51 / 1500 steps | 1.865753 [1.855085, 1.875574] | 1.866891 [1.850986, 1.872441] | 0.262738 [0.262624, 0.264812] | + +**Coefficient-only original/folded median ratios: 1.0006× GFortran, +0.9994× LFortran.** The approximately ±0.06% differences are smaller than the +observed sample spread; these measurements establish neither a meaningful +CPU benefit nor a CPU regression. The separate folded-operators/workspace +ratios were **2.9067×** and **7.1055×**, respectively, on these cases. +Do not compare absolute compiler speeds or CPU/GPU times across their +different step counts/workloads. + +Every original/folded pair had **byte-identical final fields**. Every +workspace/operator pair passed the unchanged full-field tolerances +(`rtol=1e-6`, `atol=1e-8`): maximum absolute / relative L2 differences were +`2.98e-8 / 6.64e-8` (GFortran) and `5.96e-8 / 6.90e-8` (LFortran). +All fields were finite, evolved beyond the comparator tolerances, and +preserved the initial scalar boundaries exactly. Initial setup medians +were only 87–90 µs (GFortran) and 163–175 µs (LFortran); this is the driver's +whole setup region, **not** an isolated coefficient-constructor microbenchmark. + +For an equivalent CPU build from a source snapshot, use absolute tools and +apply SDK selection to the child **after** activation: + +```bash +# From the chosen source snapshot; use a different build directory per variant. +export LFORTRAN_ROOT=/absolute/path/to/clean/lfortran +export LFORTRAN_COMPILER="$LFORTRAN_ROOT/build/llvm11/src/bin/lfortran" +export FORMAL_BENCHMARK_ROOT=/absolute/path/to/formal +unset FPM_FC FPM_FFLAGS FPM_CC FPM_CFLAGS FPM_CXX FPM_CXXFLAGS FPM_AR FPM_LDFLAGS +pixi install --locked --manifest-path "$FORMAL_BENCHMARK_ROOT/benchmarks/pixi.toml" +fpm="$FORMAL_BENCHMARK_ROOT/benchmarks/.pixi/envs/default/bin/fpm" +cpu_build="$FORMAL_BENCHMARK_ROOT/benchmarks/2d-advection-metal/build/cpu-folded-lfortran" +mkdir -p "$cpu_build/scratch" +export TMPDIR="$cpu_build/scratch" OMP_NUM_THREADS=4 + +pixi run --locked --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ + /usr/bin/env PATH="/usr/bin:/bin:$LFORTRAN_ROOT/build/llvm11/src/bin:$LFORTRAN_ROOT/.pixi/envs/llvm11/bin" \ + LFORTRAN_CC=/usr/bin/clang LFORTRAN_LINKER_PATH=/usr/bin \ + "$fpm" run --compiler "$LFORTRAN_COMPILER" --c-compiler /usr/bin/clang \ + --build-dir "$cpu_build" \ + --flag "--cpp --realloc-lhs-arrays --separate-compilation --fast" \ + --example 2D-advection-diffusion --runner /bin/cp -- "$cpu_build/formal.cpu" \ + > "$cpu_build/build.log" 2>&1 + +pixi run --locked --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ + /usr/bin/env OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 VECLIB_MAXIMUM_THREADS=1 \ + "$cpu_build/formal.cpu" --solver operators --order 4 --cells 51,51 \ + --steps 1500 --dt 1e-4 --no-csv --dump-input "$cpu_build/input.fadv2d" \ + --dump-final-raw "$cpu_build/final.fadv2d" > "$cpu_build/run.log" 2>&1 +``` + +For GFortran use the absolute +`$LFORTRAN_ROOT/.pixi/envs/llvm11/bin/gfortran`, its release flags above, and +2000 steps. Use the same flags and frozen inputs for the coefficient control. +Repeat with `--solver workspace` in the folded executable, interleave the +three variants and compare the complete exported fields. Neither building +nor a single wall-clock sample alone reproduces the controlled comparison. + +### Clean-main validation status + +The managed native CPU wrapper, missing-permission/dependency/mode +preflights, raw comparator's 54 checks in both normal Python and `-O`, CSV +self-test, and the clean LFortran `--fast` FPM suite (**32/32, zero skipped**) +passed. The actual Formal wrapper compiled and linked a Metal executable +with the clean compiler and real SDK. Its real guard's invalid-timeout +preflight was tested **before lock acquisition**, without GPU execution; +Metal/Foundation linkage and embedded kernels were separately verified. + +Shared-input CPU checks on nonsquare Gaussian 17×13/9 and nonzero, varying +boundary 13×11/3 cases passed against the unchanged native CPU oracle with +both managed Fortran compilers, including meaningful evolution and exact +boundaries. These are **CPU-only** checks. Clean-main GPU correctness, +completed GPU telemetry, and fresh 2048²/1000 timings remain pending guarded +execution; the historical table below must not be relabeled as fresh data. + +Native currently reports wall time and completed-command-buffer `gpu_busy`, +not dynamic dispatch/copy counters. Its unchanged source schedules eight +dispatches per step, eight steps per commit, and at most two command buffers +in flight. Any such native counts must be labeled **source-derived**, not +measured telemetry. Formal's runtime reports actual dispatch, commit, +completion, synchronization, copy and allocation counters. + +### Performance-portability interpretation + +The coefficient change evaluates `Ap(i,j) = -A(m+1-i,n+1-j)` in one +assignment instead of two passes. It preserves values, allocation shape and +Fortran storage order; it is not a transpose or physical layout switch. +Its CPU cost is separate from workspace allocation reuse and GPU launch +reduction. Likewise, LFortran #15506 maps adjacent GPU threads to the +contiguous subscript without requiring a Fortran array-layout change. + +Unmeasured cases include other CPUs/compilers, larger sustained CPU grids, +isolated constructor costs, and clean-main GPU/native parity. A future +automatic representation choice would need to preserve storage association, +array element order, C interoperability and observable layout (`SEQUENCE`, +`BIND(C)`, `EQUIVALENCE`, `TRANSFER`, C pointers). No such extension, +conditional CPU implementation or speculative CPU fix is included. +Bandwidth percentages from an assumed bytes-per-cell model, if discussed, +are not DRAM-counter measurements or evidence of hardware occupancy. + ## Historical paired results The original measurement used three interleaved pairs with the configuration @@ -132,6 +294,11 @@ Those values are historical verification samples, not a promise of a fixed runtime across machines. The native result was within the measured Formal throughput, with largest elementwise difference `5.96046448e-8` and relative L2 difference `5.98075916e-9`. +The compiler was snapshot `3d4c7941`, which included the now-merged GPU work +and an **unmerged driver-selection patch**. In particular, 2.195701 s is not +a fresh result from clean merged main. The earlier coefficient-only GPU +51²/500 comparison, 9.603519 → 4.057951 s (2.3666×), is also historical; +it was not automatically rerun for this CPU study. ## Index experiment diff --git a/benchmarks/2d-advection-metal/README.md b/benchmarks/2d-advection-metal/README.md index 59aa9c4..62804bb 100644 --- a/benchmarks/2d-advection-metal/README.md +++ b/benchmarks/2d-advection-metal/README.md @@ -19,8 +19,17 @@ accumulator is updated in the same `k1 + 2*k2 + 2*k3 + k4` order. ## Build and run -The wrappers use the LFortran `llvm11` Pixi environment. They accept an -explicit `LFORTRAN_ROOT`; otherwise they look for a sibling checkout at +The optional benchmark tooling uses two managed environments: + +- an already-built LFortran checkout's `llvm11` environment for its compiler, + Python and GFortran; +- [`../pixi.toml`](../pixi.toml) and [`../pixi.lock`](../pixi.lock) for + **FPM 0.13.0** and its runtime dependencies. Clean LFortran main does not + declare FPM. This small, Apple-Silicon-only benchmark environment does not + change ordinary Formal/FPM library use or the supported-compiler CI matrix. + +Install Pixi and the real Xcode Command Line Tools first. The wrappers accept +an explicit `LFORTRAN_ROOT`; otherwise they look for a sibling checkout at `../lfortran` and fail with an actionable message if it is absent: ```bash @@ -28,10 +37,17 @@ LFORTRAN_ROOT=/path/to/lfortran \ bash benchmarks/2d-advection-metal/run.sh ``` -The wrapper re-enters Pixi with -`pixi run --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11`, selects the -configured macOS SDK compiler with `xcrun`, and builds -`build/formal_advection`. With no arguments it runs a small CPU smoke case. +The wrappers re-enter +`pixi run --locked --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11`. +Their shared [`../toolchain.sh`](../toolchain.sh) selects the **real SDK +compiler after activation**, through a child-process `PATH`, `CC`, and +linker settings. This matters on clean LLVM 11 main: its Metal runtime invokes +`clang` through `PATH`; setting `CC` alone leaves Pixi's old Clang compiling +new Foundation headers and fails. There is no fake clang wrapper, compiler +patch, or change to the LFortran manifest. + +The native wrapper builds `build/formal_advection`; it does not need FPM. +With no arguments it runs a small CPU smoke case. The native executable reads `formal_advection.metal` from its current working directory; `run.sh` changes to this directory before invoking it, so it cannot accidentally load a stale shader. @@ -41,7 +57,7 @@ Metal is opt-in and doubly guarded: ```bash FA_ALLOW_METAL=1 LFORTRAN_ROOT=/path/to/lfortran \ bash benchmarks/2d-advection-metal/run.sh \ - --mode metal --cells 257,257 --steps 30000 + --mode metal --cells 17,13 --steps 9 --index-mode 2d --tg 128,1 ``` Metal runs invoke the existing LFortran @@ -67,30 +83,49 @@ That wrapper rejects missing permission or an unsupported platform before compilation. It defaults to 2048 x 2048 cells and 1000 steps, compiles `example/2D-advection-diffusion.F90` with LFortran, checks Metal linkage, and uses the same guarded launch path. It does not add native sources to FPM. +It installs the locked benchmark FPM environment automatically (or run +`pixi install --locked --manifest-path benchmarks/pixi.toml` explicitly). +FPM and Fortran are invoked by absolute managed paths, so selecting the SDK +does not accidentally select host FPM/Fortran copies. The default compiler is +`$LFORTRAN_ROOT/build/llvm11/src/bin/lfortran`; an explicit +`LFORTRAN_COMPILER` must be an absolute path. The build log records its version, +SHA-256, FPM path/version and SDK compiler. Julienne 4.1.1 is still fetched +through the unchanged `fpm.toml`. See [`PERFORMANCE.md`](./PERFORMANCE.md#reproduction) for the paired shared-input recipe and the explicit native launch geometry used in the historical comparison. +Clean-main compile/link and CPU checks have now passed with this tooling. +GPU runtime validation and timings from that clean compiler are still +pending; a successful Metal link is not a GPU execution result. See the +[controlled CPU measurements](./PERFORMANCE.md#controlled-cpu-measurements) +for the measured, workload-specific coefficient and workspace effects. + ## Known Fortran compiler limitations The standalone C/Metal reference is independent of the Fortran compiler. -The Fortran example has two known compiler-version failures in -[branch CI](https://github.com/certik/formal/actions/runs/39102014517): +The Fortran example has known compiler-version failures in +[the combined draft's initial CI](https://github.com/BerkeleyLab/formal/actions/runs/39135296334): - CI's Flang 21 crashes while lowering `do concurrent ... local(rhs, total)`. -- LFortran 0.64 fails LLVM type verification for `numerical_arguments_t` +- LFortran 0.64–0.66 fail LLVM type verification for `numerical_arguments_t` while compiling the example's argument and raw-input helpers. -Both failures occur before the solver runs. Current LFortran builds and the -other completed Flang configurations compile the portable example; real GPU -comparisons verify evolved full fields and unchanged boundaries. These older +Both failures occur before the solver runs. Current Linux LFortran builds and +the other completed Flang configurations compile the portable example; +historical real GPU comparisons verified evolved full fields and unchanged +boundaries. Fresh clean-main GPU execution remains pending. These older compiler limitations are recorded without changing the numerical method or disabling their CI checks. -That CI run also failed to build the Z3 dependency on macOS 14. The two Intel -macOS jobs were cancelled during LLVM installation, without a more specific -reported reason. None of these three jobs reached Formal compilation or tests; -they do not provide validation of those platforms. +The initial combined-draft CI passes 11 of 17 compiler checks, plus its +metadata check. Its macOS 15/26 LFortran 0.67 bottles require +`libz3.5.1.dylib`, but the installed Z3 is 5.3.1: even `lfortran --version` +fails, before Formal compilation. Linux LFortran 0.67/current and the other +Flang checks pass. These packaging failures, the older branch's macOS 14 Z3 +build failure, and its Intel macOS jobs cancelled during LLVM installation +do not validate those platforms. All remain document-only limitations: +no version drops, workflow repairs, test skips, or numerical workarounds. ## Provenance and licensing diff --git a/benchmarks/2d-advection-metal/run.sh b/benchmarks/2d-advection-metal/run.sh index b170abf..612f3f5 100755 --- a/benchmarks/2d-advection-metal/run.sh +++ b/benchmarks/2d-advection-metal/run.sh @@ -2,77 +2,61 @@ set -euo pipefail here="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" -formal_root="$(cd "$here/../.." && pwd)" -lfortran_root="${LFORTRAN_ROOT:-"$formal_root/../lfortran"}" - -if [[ ! -f "$lfortran_root/pixi.toml" ]]; then - echo "LFortran checkout not found at $lfortran_root; set LFORTRAN_ROOT to its absolute path." >&2 - exit 2 -fi -lfortran_root="$(cd "$lfortran_root" && pwd)" -if [[ "${PIXI_PROJECT_ROOT:-}" != "$lfortran_root" || - "${PIXI_ENVIRONMENT_NAME:-}" != llvm11 ]]; then - exec pixi run --manifest-path "$lfortran_root/pixi.toml" -e llvm11 \ - bash "$here/run.sh" "$@" -fi - if [[ "$(uname -s)" != Darwin || "$(uname -m)" != arm64 ]]; then echo "the native Metal benchmark requires an Apple Silicon Mac" >&2 exit 1 fi -gpu_guard="$lfortran_root/examples/nbody/gpu_run.sh" -if [[ ! -f "$gpu_guard" ]]; then - echo "LFortran GPU guard not found at $gpu_guard; use a complete LFortran checkout." >&2 - exit 2 -fi - -build_dir="$here/build" -binary="$build_dir/formal_advection" -mkdir -p "$build_dir" - -if [[ -n "${CC:-}" ]]; then - sdk_clang="$CC" -else - sdk_clang="$(xcrun --sdk macosx --find clang)" -fi -sdkroot="${SDKROOT:-$(xcrun --sdk macosx --show-sdk-path)}" -"$sdk_clang" -Wall -Wextra -Werror -O2 -fobjc-arc -isysroot "$sdkroot" \ - -framework Foundation -framework Metal \ - "$here/formal_advection.m" -o "$binary" - if [[ $# -eq 0 ]]; then set -- --mode cpu --cells 11,11 --steps 2 --dt 1e-4 fi mode=cpu mode_seen=0 -for ((i = 1; i <= $#; ++i)); do - if [[ "${!i}" == "--mode" ]]; then - if ((i == $#)); then - echo "--mode must be followed by cpu or metal" >&2 - exit 2 - fi +expect_mode=0 +for arg in "$@"; do + if ((expect_mode)); then + mode="$arg" + expect_mode=0 + elif [[ "$arg" == "--mode" ]]; then if ((mode_seen)); then echo "--mode may be specified only once" >&2 exit 2 fi - j=$((i + 1)) - mode="${!j}" + expect_mode=1 mode_seen=1 fi done +if ((expect_mode)); then + echo "--mode must be followed by cpu or metal" >&2 + exit 2 +fi if [[ "$mode" != cpu && "$mode" != metal ]]; then echo "--mode must be cpu or metal" >&2 exit 2 fi +if [[ "$mode" == metal && "${FA_ALLOW_METAL:-0}" != 1 ]]; then + echo "Metal mode is guarded; set FA_ALLOW_METAL=1 to proceed." >&2 + exit 1 +fi + +source "$here/../toolchain.sh" +benchmark_toolchain "$here/run.sh" "$@" +gpu_guard="$lfortran_root/examples/nbody/gpu_run.sh" +if [[ ! -f "$gpu_guard" ]]; then + echo "LFortran GPU guard not found at $gpu_guard; use a complete LFortran checkout." >&2 + exit 2 +fi + +build_dir="$here/build" +binary="$build_dir/formal_advection" +mkdir -p "$build_dir" +"$sdk_clang" -Wall -Wextra -Werror -O2 -fobjc-arc -isysroot "$sdkroot" \ + -framework Foundation -framework Metal \ + "$here/formal_advection.m" -o "$binary" cd "$here" if [[ "$mode" == metal ]]; then - if [[ "${FA_ALLOW_METAL:-0}" != 1 ]]; then - echo "Metal mode is guarded; set FA_ALLOW_METAL=1 to proceed." >&2 - exit 1 - fi - exec bash "$gpu_guard" "$binary" "$@" + exec /bin/bash "$gpu_guard" "$binary" "$@" fi exec "$binary" "$@" diff --git a/benchmarks/2d-advection-metal/run_formal.sh b/benchmarks/2d-advection-metal/run_formal.sh index b11b375..8f373fa 100755 --- a/benchmarks/2d-advection-metal/run_formal.sh +++ b/benchmarks/2d-advection-metal/run_formal.sh @@ -2,20 +2,6 @@ set -euo pipefail here="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" -formal_root="$(cd "$here/../.." && pwd)" -lfortran_root="${LFORTRAN_ROOT:-"$formal_root/../lfortran"}" - -if [[ ! -f "$lfortran_root/pixi.toml" ]]; then - echo "LFortran checkout not found at $lfortran_root; set LFORTRAN_ROOT to its absolute path." >&2 - exit 2 -fi -lfortran_root="$(cd "$lfortran_root" && pwd)" -if [[ "${PIXI_PROJECT_ROOT:-}" != "$lfortran_root" || - "${PIXI_ENVIRONMENT_NAME:-}" != llvm11 ]]; then - exec pixi run --manifest-path "$lfortran_root/pixi.toml" -e llvm11 \ - bash "$here/run_formal.sh" "$@" -fi - if [[ "$(uname -s)" != Darwin || "$(uname -m)" != arm64 ]]; then echo "the Formal Metal benchmark requires an Apple Silicon Mac" >&2 exit 1 @@ -25,8 +11,13 @@ if [[ "${FA_ALLOW_METAL:-0}" != 1 ]]; then exit 1 fi -compiler="${LFORTRAN_COMPILER:-${LFORTRAN_BUILD_DIR:-"$lfortran_root/build/llvm11"}/src/bin/lfortran}" +source "$here/../toolchain.sh" +benchmark_toolchain "$here/run_formal.sh" "$@" gpu_guard="$lfortran_root/examples/nbody/gpu_run.sh" +if [[ "$compiler" != /* ]]; then + echo "LFORTRAN_COMPILER must be an absolute path." >&2 + exit 2 +fi if [[ ! -x "$compiler" ]]; then echo "LFortran compiler not found at $compiler; build it with the llvm11 Pixi environment." >&2 exit 2 @@ -36,16 +27,22 @@ if [[ ! -f "$formal_root/fpm.toml" || ! -f "$gpu_guard" ]]; then exit 2 fi +benchmark_require_fpm compiler_hash="$(shasum -a 256 "$compiler" | awk '{print $1}')" build_dir="$here/build/formal/compiler-${compiler_hash:0:16}" mkdir -p "$build_dir" +printf 'compiler=%s\ncompiler_sha256=%s\nfpm=%s\nSDK clang=%s\n' \ + "$compiler" "$compiler_hash" "$benchmark_fpm" "$sdk_clang" +"$compiler" --version +"$benchmark_fpm" --version cd "$formal_root" -fpm run --compiler "$compiler" --build-dir "$build_dir" \ +"$benchmark_fpm" run --compiler "$compiler" --build-dir "$build_dir" \ + --c-compiler "$sdk_clang" \ --flag "--cpp --realloc-lhs-arrays --separate-compilation --gpu=metal" \ --c-flag "-DFORMAL_ENABLE_GPU_BENCHMARK -I$lfortran_root/src/libasr/runtime" \ - --example 2D-advection-diffusion --runner cp -- "$build_dir/formal.metal" -if ! otool -L "$build_dir/formal.metal" | grep -q '/Metal.framework/'; then + --example 2D-advection-diffusion --runner /bin/cp -- "$build_dir/formal.metal" +if ! /usr/bin/otool -L "$build_dir/formal.metal" | grep -q '/Metal.framework/'; then echo "FPM did not produce a Metal-linked executable." >&2 exit 1 fi @@ -55,5 +52,5 @@ if [[ $# -eq 0 ]]; then fi export LFORTRAN_GPU_PROFILE=1 GPU_TIMEOUT="${GPU_TIMEOUT:-30}" \ - bash "$gpu_guard" "$build_dir/formal.metal" --gpu-profile "$@" + /bin/bash "$gpu_guard" "$build_dir/formal.metal" --gpu-profile "$@" printf 'Artifacts: %s\n' "$build_dir" diff --git a/benchmarks/pixi.lock b/benchmarks/pixi.lock new file mode 100644 index 0000000..8081857 --- /dev/null +++ b/benchmarks/pixi.lock @@ -0,0 +1,99 @@ +version: 7 +platforms: +- name: osx-arm64 + virtual-packages: + - __unix=0=0 + - __osx=13.0 + - __archspec=0=m1 +environments: + default: + channels: + - url: https://conda.anaconda.org/conda-forge/ + packages: + osx-arm64: + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/_openmp_mutex-4.5-8_kmp_llvm.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/fpm-0.13.0-h20b5d57_2.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgcc-16.2.0-h3cf6597_7.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran-16.2.0-h07b0088_7.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran5-16.2.0-hdb7a957_7.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/llvm-openmp-23.1.3-hdb3d66b_1.conda +packages: +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/_openmp_mutex-4.5-8_kmp_llvm.conda + build_number: 8 + sha256: c7c98ce74f6a7ff25aaa4706d06fa41d63a333add4d697b73aa4d8d2d7ad7651 + md5: a2d706b4a7d603524133037c34f7fb76 + depends: + - llvm-openmp >=9.0.1 + license: BSD-3-Clause + license_family: BSD + run_exports: + weak: + - _openmp_mutex >=4.5 + size: 8016 + timestamp: 1788046437162 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/fpm-0.13.0-h20b5d57_2.conda + sha256: ccbe33259643ceebaabacbb28487fa9daaca52e1d0cf52af90ae4dc88dc78fe6 + md5: 642d5ecc03f178a86d61358c2c98ec27 + depends: + - __osx >=11.0 + - libgcc >=15 + - libgfortran + - libgfortran5 >=15.3.0 + - llvm-openmp >=21.1.8 + license: MIT + license_family: MIT + run_exports: {} + size: 883579 + timestamp: 1791010730175 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgcc-16.2.0-h3cf6597_7.conda + sha256: deacf62bd6cbbe05eb085dd69a2ca120053597c37066d53c8f7cf6815f26e9eb + md5: f100c8baeb0c816052b0d43948555bf6 + depends: + - _openmp_mutex + constrains: + - libgcc-ng ==16.2.0=*_7 + - libgomp 16.2.0 7 + license: GPL-3.0-only WITH GCC-exception-3.1 + license_family: GPL + run_exports: {} + size: 365973 + timestamp: 1790301995634 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran-16.2.0-h07b0088_7.conda + sha256: a5b4b63d0ef4c7c59c37225b9f0307b5cd2a3a8499b2a40b6ddd2b9f072712a8 + md5: 2efa9518e126c501a7c7d638b2fa3a21 + depends: + - libgfortran5 16.2.0 hdb7a957_7 + constrains: + - libgfortran-ng ==16.2.0=*_7 + license: GPL-3.0-only WITH GCC-exception-3.1 + license_family: GPL + run_exports: {} + size: 99745 + timestamp: 1790302083840 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libgfortran5-16.2.0-hdb7a957_7.conda + sha256: 19297e5cf39de976cdd1a46499a4e15fdf54d249d3231f3b2a9d43af3be0bc66 + md5: c9d80dd8775b8fe76e825d63d8ccaf1b + depends: + - libgcc >=16.2.0 + constrains: + - libgfortran 16.2.0 + license: GPL-3.0-only WITH GCC-exception-3.1 + license_family: GPL + run_exports: {} + size: 554793 + timestamp: 1790302001166 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/llvm-openmp-23.1.3-hdb3d66b_1.conda + sha256: eec66f3d0a0c428411e22a59c4726069ccef81a46864876b6e270c8c254d4045 + md5: 5b7efd98534e6722de8d223b8654185b + depends: + - __osx >=11.0 + constrains: + - intel-openmp <0.0a0 + - openmp 23.1.3|23.1.3.* + license: Apache-2.0 WITH LLVM-exception + license_family: APACHE + run_exports: + strong: + - llvm-openmp >=23.1.3 + size: 287944 + timestamp: 1791950643503 diff --git a/benchmarks/pixi.toml b/benchmarks/pixi.toml new file mode 100644 index 0000000..65f171a --- /dev/null +++ b/benchmarks/pixi.toml @@ -0,0 +1,8 @@ +[workspace] +name = "formal-benchmarks" +description = "FPM for the optional Apple Silicon advection benchmarks" +channels = ["conda-forge"] +platforms = ["osx-arm64"] + +[dependencies] +fpm = "==0.13.0" diff --git a/benchmarks/toolchain.sh b/benchmarks/toolchain.sh new file mode 100644 index 0000000..568ac96 --- /dev/null +++ b/benchmarks/toolchain.sh @@ -0,0 +1,51 @@ +#!/usr/bin/env bash + +# Source this helper, then call benchmark_toolchain with the entry script and arguments. +benchmark_toolchain() { + local entry_script="$1" + shift + formal_root="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" + lfortran_root="${LFORTRAN_ROOT:-"$formal_root/../lfortran"}" + if [[ ! -f "$lfortran_root/pixi.toml" ]]; then + echo "LFortran checkout not found at $lfortran_root; set LFORTRAN_ROOT to its absolute path." >&2 + exit 2 + fi + lfortran_root="$(cd "$lfortran_root" && pwd)" + benchmark_pixi="$(command -v pixi)" || { + echo "Pixi is required for the benchmark tools; install Pixi before running this script." >&2 + exit 2 + } + if [[ "${PIXI_PROJECT_ROOT:-}" != "$lfortran_root" || + "${PIXI_ENVIRONMENT_NAME:-}" != llvm11 ]]; then + export LFORTRAN_ROOT="$lfortran_root" + exec "$benchmark_pixi" run --locked --manifest-path "$lfortran_root/pixi.toml" -e llvm11 \ + /bin/bash "$entry_script" "$@" + fi + + compiler="${LFORTRAN_COMPILER:-"$lfortran_root/build/llvm11/src/bin/lfortran"}" + benchmark_fpm="$formal_root/benchmarks/.pixi/envs/default/bin/fpm" + benchmark_python="$lfortran_root/.pixi/envs/llvm11/bin/python" + benchmark_gfortran="$lfortran_root/.pixi/envs/llvm11/bin/gfortran" + + # LFortran's Metal runtime invokes clang through PATH, not CC alone. + # Select the actual SDK tools after Pixi activation, without replacing + # the explicitly selected Fortran, FPM or Python executables. + if ! sdk_clang="$(/usr/bin/xcrun --sdk macosx --find clang)" || + ! sdk_clangxx="$(/usr/bin/xcrun --sdk macosx --find clang++)" || + ! sdkroot="$(/usr/bin/xcrun --sdk macosx --show-sdk-path)" || + [[ ! -x "$sdk_clang" || ! -d "$sdkroot" ]]; then + echo "The macOS SDK is required; install Xcode Command Line Tools." >&2 + exit 2 + fi + export PATH="$(dirname "$sdk_clang"):/usr/bin:/bin:$lfortran_root/build/llvm11/src/bin:$lfortran_root/.pixi/envs/llvm11/bin:$PATH" + export CC="$sdk_clang" CXX="$sdk_clangxx" LFORTRAN_CC="$sdk_clang" + export LFORTRAN_LINKER_PATH=/usr/bin SDKROOT="$sdkroot" +} + +benchmark_require_fpm() { + "$benchmark_pixi" install --locked --manifest-path "$formal_root/benchmarks/pixi.toml" + if [[ ! -x "$benchmark_fpm" ]]; then + echo "Managed FPM is missing; run pixi install --locked --manifest-path benchmarks/pixi.toml." >&2 + exit 2 + fi +} From b8d5c219038c14f9f6d474c4b7185e9037fe80dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sun, 18 Oct 2026 09:43:44 -0600 Subject: [PATCH 15/18] docs: record clean-main GPU performance measurements --- benchmarks/2d-advection-metal/PERFORMANCE.md | 118 +++++++++++++++++-- benchmarks/2d-advection-metal/README.md | 15 ++- 2 files changed, 114 insertions(+), 19 deletions(-) diff --git a/benchmarks/2d-advection-metal/PERFORMANCE.md b/benchmarks/2d-advection-metal/PERFORMANCE.md index cec4a69..edf3a54 100644 --- a/benchmarks/2d-advection-metal/PERFORMANCE.md +++ b/benchmarks/2d-advection-metal/PERFORMANCE.md @@ -187,6 +187,9 @@ allocation/copy is still inside the solve. observed sample spread; these measurements establish neither a meaningful CPU benefit nor a CPU regression. The separate folded-operators/workspace ratios were **2.9067×** and **7.1055×**, respectively, on these cases. +Coefficient folding is CPU-neutral within this sample spread; the workspace +ratios compare workspace against **folded operators**, not old versus folded +coefficients. Do not compare absolute compiler speeds or CPU/GPU times across their different step counts/workloads. @@ -237,31 +240,118 @@ Repeat with `--solver workspace` in the folded executable, interleave the three variants and compare the complete exported fields. Neither building nor a single wall-clock sample alone reproduces the controlled comparison. -### Clean-main validation status +### Tooling and CPU validation The managed native CPU wrapper, missing-permission/dependency/mode preflights, raw comparator's 54 checks in both normal Python and `-O`, CSV self-test, and the clean LFortran `--fast` FPM suite (**32/32, zero skipped**) passed. The actual Formal wrapper compiled and linked a Metal executable -with the clean compiler and real SDK. Its real guard's invalid-timeout -preflight was tested **before lock acquisition**, without GPU execution; -Metal/Foundation linkage and embedded kernels were separately verified. +with the clean compiler and real SDK. Before runtime validation, its real +guard's invalid-timeout preflight was tested **before lock acquisition**, +without GPU execution; Metal/Foundation linkage and embedded kernels were +separately verified. The subsequent guarded executions are reported below. Shared-input CPU checks on nonsquare Gaussian 17×13/9 and nonzero, varying boundary 13×11/3 cases passed against the unchanged native CPU oracle with both managed Fortran compilers, including meaningful evolution and exact -boundaries. These are **CPU-only** checks. Clean-main GPU correctness, -completed GPU telemetry, and fresh 2048²/1000 timings remain pending guarded -execution; the historical table below must not be relabeled as fresh data. +boundaries. These checks were **CPU-only**; the following section records +the separate, subsequently completed GPU executions. + +## Clean-main paired GPU results + +Measured on **18 October 2026**, using Formal source/tooling head +`70e73ee51f72ebaca16544d834bdfd98a9da66a4` and the unpatched clean LFortran +`2247d35082ea190408486b93c12ec2ead4771735` / LLVM 11.1.0 compiler identified +above. The machine was the same M4 Max/macOS 15.8 system: AC power, +`powermode=0`, no core affinity or clock locking, load averages about +2.64–2.66. The coordinator ran the prepared procedure successfully; this +documentation follow-up only audited the saved artifacts, without rerunning +any solver or benchmark. + +All **ten GPU executions** passed the unchanged shared guard, with the +actual executable immediately after it, a 30 s timeout and 4096 MiB RSS +limit per executable. The two small Formal/native pairs ran first: + +- **Gaussian 17×13/9**, shared raw input: Formal/native maximum difference + `2.32830644e-10`; each GPU final agreed with the native CPU oracle to + `7.45058060e-9` maximum absolute difference. +- **Nonzero, varying boundary 13×11/3**: Formal, native GPU and CPU-oracle + finals were identical. All three evolved meaningfully, rather than + merely returning identical initial fields. + +Both small pairs reported positive completed GPU time and passed full-field +finiteness, evolution and exact fixed-boundary checks at the unchanged +`rtol=1e-6`, `atol=1e-8`. + +The substantial case was the **actual Gaussian, 2048×2048 cells, 1000 steps**, +with exact FP32 `dt=6.20126741068816e-8`, one shared raw input, and native +**direct 2D / 128×1** geometry. The first Formal run exported that input; +the later Formal runs and every native run consumed the same immutable +snapshot. Three interleaved pairs ran Formal/native, native/Formal, +Formal/native. Each executable performed an untimed one-step warmup and +restored its initial state before solve timing. + +| Implementation | Median solve wall [min, max], s | Median completed GPU busy [min, max], s | +|---|---:|---:| +| Formal + clean-main LFortran workspace | **2.167132** [2.161986, 2.167841] | **2.158323** [2.152374, 2.159060] | +| Native C + Metal | **2.188441** [2.184408, 2.203517] | **2.187208** [2.183102, 2.202095] | + +Native/Formal median ratios are **1.00983× wall** and **1.01338× GPU busy**. +Both implementations therefore solve this case in about 2.2 s; this is +workload/toolchain-specific evidence, not a universal speed or occupancy +limit. Timings exclude compilation, initialization, warmup/JIT and field I/O. +All six large solves stayed in the specified 1–5 s window. + +Every large pair passed the complete field comparison: maximum absolute +difference **`5.96046448e-8`**, relative L2 difference +**`5.98075916e-9`**, finite values and exactly preserved initial boundaries. +Both implementations changed the scalar by about `4.17232513e-4` in max norm +and `2.8552e-4` in relative L2, well beyond the evolution thresholds. +Within each implementation, final files were byte-identical across all +three repetitions; Formal and native finals were not byte-identical to +each other. The large case was compared between GPU implementations; +CPU-oracle validation was performed on the two bounded small cases. + +Formal reported the following **observed runtime counters in each large +solve**, not just the expected `8 * steps` stage-kernel count: + +| Counter | Per solve | +|---|---:| +| Dispatches | 8,003 | +| Commits / completed command buffers | 126 / 126 | +| Host syncs / waits | 2 / 124 | +| Host-to-device / device-to-host memcpy bytes | 201,654,400 / 117,637,200 | +| Buffer allocations / allocated bytes | 5 / 84,017,200 | +| JIT time inside solve | 0 s | + +These byte counts describe the runtime's logical copies, **not physical +DRAM traffic measurements**. Five runtime allocations also preclude an +allocation-free claim. Native currently reports wall time and completed-command-buffer `gpu_busy`, not dynamic dispatch/copy counters. Its unchanged source schedules eight dispatches per step, eight steps per commit, and at most two command buffers -in flight. Any such native counts must be labeled **source-derived**, not -measured telemetry. Formal's runtime reports actual dispatch, commit, -completion, synchronization, copy and allocation counters. +in flight: 8,000 dispatches, 125 commits/waits and a 16,810,000-byte final +shared-buffer-to-host copy for this case. These native counts are +**source-derived**, not measured telemetry. + +Identities of the measured artifacts (not guarantees of bit-identical builds +on other machines): + +```text +Formal Metal binary SHA256: +3b0f83ac19fb8858cb59863dffd5da2c882ccdcb801430f6007b008943485668 +Native binary SHA256: +7873f33da20935ab50dae00c5a319d0f43301d1f837d54d37804b9eca2de485d +Shared 2048x2048/1000 raw input SHA256: +ad304f4fce24f392965f4bca34f6d74f8a57570da185d2547537db91cdefbfc7 +``` + +Source/binary/guard/input identities and every reported sample were audited +against the saved execution and comparison logs. No historical binary or +old result was substituted. -### Performance-portability interpretation +## Performance-portability interpretation The coefficient change evaluates `Ap(i,j) = -A(m+1-i,n+1-j)` in one assignment instead of two passes. It preserves values, allocation shape and @@ -271,7 +361,7 @@ reduction. Likewise, LFortran #15506 maps adjacent GPU threads to the contiguous subscript without requiring a Fortran array-layout change. Unmeasured cases include other CPUs/compilers, larger sustained CPU grids, -isolated constructor costs, and clean-main GPU/native parity. A future +isolated constructor costs, and GPU/native parity on other workloads. A future automatic representation choice would need to preserve storage association, array element order, C interoperability and observable layout (`SEQUENCE`, `BIND(C)`, `EQUIVALENCE`, `TRANSFER`, C pointers). No such extension, @@ -298,7 +388,9 @@ The compiler was snapshot `3d4c7941`, which included the now-merged GPU work and an **unmerged driver-selection patch**. In particular, 2.195701 s is not a fresh result from clean merged main. The earlier coefficient-only GPU 51²/500 comparison, 9.603519 → 4.057951 s (2.3666×), is also historical; -it was not automatically rerun for this CPU study. +it was not rerun for this study. The fresh table above establishes clean-main +application behavior; the difference from the historical table is not a +controlled old/new compiler speedup comparison. ## Index experiment diff --git a/benchmarks/2d-advection-metal/README.md b/benchmarks/2d-advection-metal/README.md index 62804bb..2124e0c 100644 --- a/benchmarks/2d-advection-metal/README.md +++ b/benchmarks/2d-advection-metal/README.md @@ -95,9 +95,12 @@ See [`PERFORMANCE.md`](./PERFORMANCE.md#reproduction) for the paired shared-input recipe and the explicit native launch geometry used in the historical comparison. -Clean-main compile/link and CPU checks have now passed with this tooling. -GPU runtime validation and timings from that clean compiler are still -pending; a successful Metal link is not a GPU execution result. See the +Clean-main compile/link, CPU checks and ten guarded GPU executions have now +passed with this tooling. The substantial 2048²/1000 Gaussian comparison +measured median solve walls of **2.167132 s Formal / 2.188441 s native**, +with full-field, evolution and fixed-boundary validation. These are fresh +clean-main results, not merely a successful Metal link. See the +[paired GPU measurements](./PERFORMANCE.md#clean-main-paired-gpu-results) and [controlled CPU measurements](./PERFORMANCE.md#controlled-cpu-measurements) for the measured, workload-specific coefficient and workspace effects. @@ -112,9 +115,9 @@ The Fortran example has known compiler-version failures in while compiling the example's argument and raw-input helpers. Both failures occur before the solver runs. Current Linux LFortran builds and -the other completed Flang configurations compile the portable example; -historical real GPU comparisons verified evolved full fields and unchanged -boundaries. Fresh clean-main GPU execution remains pending. These older +the other completed Flang configurations compile the portable example. +Historical and fresh clean-main GPU comparisons verify evolved full fields +and unchanged boundaries. These older compiler limitations are recorded without changing the numerical method or disabling their CI checks. From a5a3db72c1aa21a76c7cab9e0517fd2eef67752e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sun, 18 Oct 2026 11:20:37 -0600 Subject: [PATCH 16/18] Check advection workspace bounds and finite timesteps --- example/2D-advection-diffusion.F90 | 1 + example/formal_advection_workspace.F90 | 41 ++++++++++++++++------ example/validate-advection-inputs.f90 | 47 ++++++++++++++++++++++++++ 3 files changed, 78 insertions(+), 11 deletions(-) create mode 100644 example/validate-advection-inputs.f90 diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90 index 79cef80..d35803e 100644 --- a/example/2D-advection-diffusion.F90 +++ b/example/2D-advection-diffusion.F90 @@ -398,6 +398,7 @@ subroutine validate_numerical_arguments(numerical_arguments) if (any(numerical_arguments%cells < 2*numerical_arguments%order + 1)) error stop "cells must satisfy cells >= 2*order + 1" if (numerical_arguments%steps <= 0) error stop "steps must be positive" if (.not. allocated(numerical_arguments%input_path)) then + call fa_validate_working_set(numerical_arguments%cells(1), numerical_arguments%cells(2)) if (.not. numerical_arguments%dt_supplied .and. numerical_arguments%dt <= 0E0) error stop "dt must be positive" if (numerical_arguments%dt <= 0E0) error stop "dt must be positive" end if diff --git a/example/formal_advection_workspace.F90 b/example/formal_advection_workspace.F90 index 0bf5c0e..2376b8e 100644 --- a/example/formal_advection_workspace.F90 +++ b/example/formal_advection_workspace.F90 @@ -89,30 +89,49 @@ end subroutine fa_die integer(int64) function fa_checked_mul(a, b, what) result(product) integer(int64), intent(in) :: a, b character(len=*), intent(in) :: what - if (a /= 0_int64 .and. b > huge(0_int64) / a) then - call fa_die(trim(what) // " size overflow") + if (a < 0_int64 .or. b < 0_int64) call fa_die(trim(what) // " size must be non-negative") + if (a > 0_int64) then + if (b > huge(0_int64) / a) call fa_die(trim(what) // " size overflow") end if product = a * b end function fa_checked_mul + integer(int64) function fa_checked_add(a, b, what) result(total) + integer(int64), intent(in) :: a, b + character(len=*), intent(in) :: what + if (a < 0_int64 .or. b < 0_int64) call fa_die(trim(what) // " size must be non-negative") + if (a > huge(0_int64) - b) call fa_die(trim(what) // " size overflow") + total = a + b + end function fa_checked_add + subroutine fa_validate_working_set(nx, ny) integer, intent(in) :: nx, ny - integer(int64) :: scalar_len, xface_len, yface_len, total_bytes - scalar_len = fa_checked_mul(int(nx + 2, int64), int(ny + 2, int64), & + integer(int64) :: nx64, ny64, scalar_len, xface_len, yface_len + integer(int64) :: scalar_bytes, xface_bytes, yface_bytes, total_bytes + if (nx < 1 .or. ny < 1) call fa_die("cells must be positive") + if (nx > huge(nx) - 2 .or. ny > huge(ny) - 2) then + call fa_die("cells including halo size overflow") + end if + nx64 = int(nx, int64) + ny64 = int(ny, int64) + scalar_len = fa_checked_mul(nx64 + 2_int64, ny64 + 2_int64, & "scalar field") - xface_len = fa_checked_mul(int(nx + 1, int64), int(ny + 2, int64), & + xface_len = fa_checked_mul(nx64 + 1_int64, ny64 + 2_int64, & "x-face field") - yface_len = fa_checked_mul(int(nx + 2, int64), int(ny + 1, int64), & + yface_len = fa_checked_mul(nx64 + 2_int64, ny64 + 1_int64, & "y-face field") - total_bytes = 3_int64 * scalar_len * 4_int64 & - + 3_int64 * xface_len * 4_int64 & - + 3_int64 * yface_len * 4_int64 + scalar_bytes = fa_checked_mul(scalar_len, 12_int64, "scalar buffers") + xface_bytes = fa_checked_mul(xface_len, 12_int64, "x-face buffers") + yface_bytes = fa_checked_mul(yface_len, 12_int64, "y-face buffers") + total_bytes = fa_checked_add(scalar_bytes, xface_bytes, "working set") + total_bytes = fa_checked_add(total_bytes, yface_bytes, "working set") if (total_bytes > 512_int64 * 1024_int64 * 1024_int64) then call fa_die("estimated working set exceeds 512 MiB") end if end subroutine fa_validate_working_set subroutine fa_validate_cfg(cfg, dt_supplied) + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite type(fa_config_t), intent(inout) :: cfg logical, intent(in) :: dt_supplied real(real32) :: scale @@ -122,7 +141,7 @@ subroutine fa_validate_cfg(cfg, dt_supplied) if (cfg%steps < 1) then call fa_die("steps must be positive") end if - if (.not. (cfg%dt > 0.0_real32)) then + if (.not. ieee_is_finite(cfg%dt) .or. .not. (cfg%dt > 0.0_real32)) then call fa_die("dt must be finite and positive") end if if (.not. dt_supplied) then @@ -231,7 +250,6 @@ subroutine fa_read_raw_input(path, fields, nx, ny, steps, dt) ny = cfg%ny steps = cfg%steps dt = cfg%dt - call fa_validate_working_set(nx, ny) call fa_allocate_fields(fields, nx, ny) read(unit, iostat=ios) fields%base if (ios /= 0) call fa_die("cannot read raw base field") @@ -433,6 +451,7 @@ end function fa_div_face_y subroutine fa_allocate_fields(fields, nx, ny) type(fa_fields_t), intent(inout) :: fields integer, intent(in) :: nx, ny + call fa_validate_working_set(nx, ny) allocate(fields%base(1:nx+2, 1:ny+2)) allocate(fields%q(1:nx+2, 1:ny+2)) allocate(fields%acc(1:nx+2, 1:ny+2)) diff --git a/example/validate-advection-inputs.f90 b/example/validate-advection-inputs.f90 new file mode 100644 index 0000000..e19a4b8 --- /dev/null +++ b/example/validate-advection-inputs.f90 @@ -0,0 +1,47 @@ +! Copyright (c) 2026, The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +program validate_advection_inputs + use formal_advection_workspace, only: fa_config_t, fa_validate_cfg, fa_validate_working_set + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_positive_inf, ieee_negative_inf + implicit none + type(fa_config_t) :: cfg + character(len=64) :: mode, argument + integer :: nx, ny, ios + + ! This host-only harness never allocates fields or calls the solver. + call get_command_argument(1, mode) + select case (trim(mode)) + case ("working-set") + if (command_argument_count() /= 3) error stop "expected working-set nx ny" + call get_command_argument(2, argument) + read(argument, *, iostat=ios) nx + if (ios /= 0) error stop "invalid nx" + call get_command_argument(3, argument) + read(argument, *, iostat=ios) ny + if (ios /= 0) error stop "invalid ny" + call fa_validate_working_set(nx, ny) + case ("timestep") + if (command_argument_count() /= 2) error stop "expected timestep dt" + call get_command_argument(2, argument) + select case (trim(argument)) + case ("nan") + cfg%dt = ieee_value(cfg%dt, ieee_quiet_nan) + case ("inf", "+inf") + cfg%dt = ieee_value(cfg%dt, ieee_positive_inf) + case ("-inf") + cfg%dt = ieee_value(cfg%dt, ieee_negative_inf) + case default + read(argument, *, iostat=ios) cfg%dt + if (ios /= 0) error stop "invalid timestep" + end select + call fa_validate_cfg(cfg, .true.) + case ("defaults") + if (command_argument_count() /= 1) error stop "expected defaults" + call fa_validate_cfg(cfg, .false.) + call fa_validate_working_set(cfg%nx, cfg%ny) + case default + error stop "expected working-set, timestep or defaults" + end select + print '(a)', "input validation passed" +end program validate_advection_inputs From 1e3319fc22ca482fb07092fec43ecb6d107164b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sun, 18 Oct 2026 11:20:37 -0600 Subject: [PATCH 17/18] Share strict field validation across advection tools --- example/scripts/compare_2d_fields.py | 29 +- example/scripts/compare_2d_raw_fields.py | 121 +++++++-- example/scripts/test_advection_validation.py | 247 ++++++++++++++++++ .../validate_2d_advection_workspace.py | 181 ++++++------- 4 files changed, 450 insertions(+), 128 deletions(-) create mode 100644 example/scripts/test_advection_validation.py diff --git a/example/scripts/compare_2d_fields.py b/example/scripts/compare_2d_fields.py index ecc8ab1..2b0572a 100755 --- a/example/scripts/compare_2d_fields.py +++ b/example/scripts/compare_2d_fields.py @@ -65,9 +65,27 @@ def load_field(path: Path): return xs, ys, field -def compare_fields(reference, observed, rtol: float, atol: float): - if len(reference) != len(observed) or len(reference[0]) != len(observed[0]): +def validate_tolerances(rtol: float, atol: float): + if not all(math.isfinite(t) and t >= 0 for t in (rtol, atol)): + raise SystemExit("--rtol and --atol must be finite and non-negative") + + +def validate_fields(reference, observed): + if not reference or not observed or not reference[0] or not observed[0]: + raise SystemExit("empty field data") + if len(reference) != len(observed): raise SystemExit("field shapes differ") + ny = len(reference[0]) + for field in (reference, observed): + if any(len(row) != ny for row in field): + raise SystemExit("field shapes differ") + if not all(math.isfinite(value) for row in field for value in row): + raise SystemExit("non-finite field data") + + +def compare_fields(reference, observed, rtol: float, atol: float): + validate_tolerances(rtol, atol) + validate_fields(reference, observed) max_abs_diff = 0.0 max_rel_diff = 0.0 @@ -93,8 +111,7 @@ def compare_fields(reference, observed, rtol: float, atol: float): def summarize_field_differences(reference, observed): - if len(reference) != len(observed) or len(reference[0]) != len(observed[0]): - raise SystemExit("field shapes differ") + validate_fields(reference, observed) max_abs_diff = 0.0 max_rel_diff = 0.0 @@ -123,6 +140,7 @@ def load_and_validate(reference_path: Path, observed_path: Path): def report_comparison(reference, observed, rtol: float, atol: float, evolution: bool): + validate_tolerances(rtol, atol) if evolution: max_abs_diff, max_rel_diff, relative_l2, changed = summarize_field_differences(reference, observed) else: @@ -215,8 +233,7 @@ def main(argv): if args.reference is None or args.candidate is None: parser.error("reference and candidate CSV paths are required unless --self-test is used") - if not math.isfinite(args.rtol) or not math.isfinite(args.atol) or args.rtol < 0.0 or args.atol < 0.0: - raise SystemExit("--rtol and --atol must be finite and non-negative") + validate_tolerances(args.rtol, args.atol) run_comparison(Path(args.reference), Path(args.candidate), args.rtol, args.atol, args.evolution) return 0 diff --git a/example/scripts/compare_2d_raw_fields.py b/example/scripts/compare_2d_raw_fields.py index abf7c9c..a163882 100644 --- a/example/scripts/compare_2d_raw_fields.py +++ b/example/scripts/compare_2d_raw_fields.py @@ -3,14 +3,89 @@ import argparse import math +import struct import sys from pathlib import Path -from compare_2d_fields import report_comparison -from validate_2d_advection_workspace import parse_raw_headered, parse_raw_headerless +from compare_2d_fields import report_comparison, validate_fields, validate_tolerances + + +def parse_col_major_payload(payload: bytes, nx: int, ny: int): + if nx < 1 or ny < 1: + raise SystemExit("raw field dimensions must be positive") + count = (nx + 2) * (ny + 2) + if len(payload) != count * 4: + raise SystemExit(f"unexpected raw payload size: expected {count * 4} bytes") + values = struct.unpack(f"<{count}f", payload) + if not all(math.isfinite(value) for value in values): + raise SystemExit("non-finite raw scalar field") + return [[values[i + j * (nx + 2)] for j in range(ny + 2)] for i in range(nx + 2)] + + +def fp32_timestep(raw, label): + try: + dt = struct.unpack("= 9 and a positive step count") + if snapshot[:3] != (*cells, steps): + raise SystemExit("initial configuration mismatch") + if dt is not None and snapshot[3] != fp32_timestep(dt, "expected dt"): + raise SystemExit("initial timestep differs from expected dt") + + +def read_matching_final(path, initial_snapshot, label): + snapshot = parse_raw_headered(path, expected_mode="raw-final") + if snapshot[:4] != initial_snapshot[:4]: + raise SystemExit(f"{label} configuration differs from shared input") + return snapshot[4] def check_boundaries(initial, final, label): + validate_fields(initial, final) nx, ny = len(initial) - 2, len(initial[0]) - 2 boundary = [(i, j) for i in range(nx + 2) for j in (0, ny + 1)] boundary += [(i, j) for i in (0, nx + 1) for j in range(1, ny + 1)] @@ -19,38 +94,38 @@ def check_boundaries(initial, final, label): raise SystemExit(f"{label}: fixed boundary changed at ({i},{j})") +def check_final_field(initial, final, label, rtol, atol, evolution=False): + validate_tolerances(rtol, atol) + validate_fields(initial, final) + if evolution: + initial_scale = max(abs(value) for row in initial for value in row) + if initial_scale == 0: + raise SystemExit("--evolution requires a nonzero initial field") + if not any(value != 0 for row in final for value in row): + raise SystemExit(f"{label}: zero final field from nonzero input") + check_boundaries(initial, final, label) + if evolution: + print(f"{label} evolution:") + max_abs, _, relative_l2, _ = report_comparison(initial, final, rtol, atol, True) + if max_abs <= atol + rtol * initial_scale or relative_l2 <= rtol: + raise SystemExit(f"{label}: evolution does not exceed comparison tolerances") + + def compare_raw_fields(initial_path, formal_path, native_path, cells, steps, rtol=1e-6, atol=1e-8, evolution=False, cpu_path=None): - if len(cells) != 2 or min(cells) < 9 or steps <= 0: - raise SystemExit("expected order-4 cells >= 9 and a positive step count") - if not all(math.isfinite(t) and t >= 0 for t in (rtol, atol)): - raise SystemExit("--rtol and --atol must be finite and non-negative") + validate_tolerances(rtol, atol) initial_snapshot = parse_raw_headered(initial_path, expected_mode="raw") - final_snapshot = parse_raw_headered(formal_path, expected_mode="raw-final") + check_initial_configuration(initial_snapshot, cells, steps) nx, ny, n, dt, initial = initial_snapshot - if (nx, ny, n) != (*cells, steps): - raise SystemExit("initial configuration mismatch") - if final_snapshot[:4] != initial_snapshot[:4]: - raise SystemExit("Formal final configuration differs from shared input") - formal = final_snapshot[4] + formal = read_matching_final(formal_path, initial_snapshot, "Formal final") native = parse_raw_headerless(native_path, nx, ny) finals = [("Formal", formal), ("native", native)] if cpu_path is not None: finals.append(("CPU oracle", parse_raw_headerless(cpu_path, nx, ny))) - initial_scale = max(abs(value) for row in initial for value in row) - if evolution and initial_scale == 0: - raise SystemExit("--evolution requires a nonzero initial field") print(f"cells={nx},{ny} steps={n} dt={dt:.16e} elements={(nx + 2) * (ny + 2)}") for label, final in finals: - if evolution and not any(value != 0 for row in final for value in row): - raise SystemExit(f"{label}: zero final field from nonzero input") - check_boundaries(initial, final, label) - if evolution: - print(f"{label} evolution:") - max_abs, _, relative_l2, _ = report_comparison(initial, final, rtol, atol, True) - if max_abs <= atol + rtol * initial_scale or relative_l2 <= rtol: - raise SystemExit(f"{label}: evolution does not exceed comparison tolerances") + check_final_field(initial, final, label, rtol, atol, evolution) print("Formal versus native:") report_comparison(formal, native, rtol, atol, False) if cpu_path is not None: diff --git a/example/scripts/test_advection_validation.py b/example/scripts/test_advection_validation.py new file mode 100644 index 0000000..a9fe355 --- /dev/null +++ b/example/scripts/test_advection_validation.py @@ -0,0 +1,247 @@ +from __future__ import annotations + +import contextlib +import io +import math +from pathlib import Path +import shutil +import struct +import subprocess +import unittest +from unittest.mock import patch + +import compare_2d_fields as csv_fields +import compare_2d_raw_fields as raw_fields +import validate_2d_advection_workspace as workspace + + +def snapshot(nx=13, ny=11, steps=3, dt="0.0001", mode="raw-final", bad_value=None): + scalar = [0.5 + 0.01 * i + 0.02 * j for j in range(ny + 2) for i in range(nx + 2)] + if mode == "raw-final": + for j in range(1, ny + 1): + for i in range(1, nx + 1): + scalar[i + j * (nx + 2)] *= 0.99 + values = scalar + if mode == "raw": + values += [0.1] * ((nx + 1) * (ny + 2) + (nx + 2) * (ny + 1)) + if bad_value is not None: + index, value = bad_value + values[index] = value + header = f"FADV2D1\n{nx} {ny} {steps} {dt}\nmode {mode}\n".encode() + return header + struct.pack(f"<{len(values)}f", *values) + + +class WorkspaceValidationTests(unittest.TestCase): + cases = 0 + + @classmethod + def setUpClass(cls): + cls.scratch = Path(__file__).resolve().parents[2] / "build/formal-gpu/validate-workspace-selftest" + cls.scratch.mkdir(parents=True, exist_ok=False) + cls.addClassCleanup(shutil.rmtree, cls.scratch) + + def check_driver(self, *, initial=None, reference=None, candidate=None, native=None, + extra=(), error=None, launched=True, write_candidate=True): + initial_path = self.scratch / "input.fadv2d" + reference_path = self.scratch / "reference.fadv2d" + initial_path.write_bytes(initial if initial is not None else snapshot(mode="raw")) + reference_path.write_bytes(reference if reference is not None else snapshot()) + argv = ["validate", "--executable", str(self.scratch / "not-a-solver"), + "--gpu-runner", "/usr/bin/env", "--workdir", str(self.scratch), + "--raw-input", str(initial_path), "--reference-raw", str(reference_path)] + if native is not None: + native_path = self.scratch / "native.bin" + native_path.write_bytes(native) + argv += ["--native-final-raw", str(native_path)] + + def emit(command, cwd): + if write_candidate: + path = Path(command[command.index("--dump-final-raw") + 1]) + path.write_bytes(candidate if candidate is not None else snapshot()) + + with patch.object(workspace, "run", side_effect=emit) as launch, \ + contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(io.StringIO()): + if error is None: + self.assertEqual(workspace.main(argv + list(extra)), 0) + else: + with self.assertRaisesRegex((SystemExit, OSError), error): + workspace.main(argv + list(extra)) + self.assertEqual(launch.call_count, int(launched)) + type(self).cases += 1 + + def test_valid_metadata_and_fp32_rounding(self): + self.check_driver(extra=("--evolution",)) + self.check_driver(extra=("--cells", "13,11", "--steps", "3", "--dt", "0.0001")) + rounded = f"{raw_fields.fp32_timestep('0.0001', 'test'):.16e}" + self.check_driver(reference=snapshot(dt=rounded), candidate=snapshot(dt=rounded), + extra=("--dt", rounded)) + self.check_driver(native=snapshot().split(b"\n", 3)[3]) + self.check_driver(extra=("--rtol", "0", "--atol", "0")) + + def test_wrong_candidate_and_reference_metadata(self): + for field in ("reference", "candidate"): + for key, value in (("nx", 14), ("ny", 12), ("steps", 999), ("dt", "0.123")): + with self.subTest(field=field, key=key): + self.check_driver(**{field: snapshot(**{key: value})}, + error="configuration differs", launched=field == "candidate") + + def test_wrong_expected_metadata(self): + for extra in (("--cells", "14,11"), ("--cells", "13,12"), ("--steps", "4")): + self.check_driver(extra=extra, error="initial configuration mismatch", launched=False) + self.check_driver(extra=("--dt", "0.0002"), error="initial timestep differs", launched=False) + + def test_nonfinite_or_negative_tolerances(self): + for option in ("--rtol", "--atol"): + for value in ("nan", "inf", "-inf", "-1"): + with self.subTest(option=option, value=value): + self.check_driver(extra=(f"{option}={value}",), + error="finite and non-negative", launched=False) + with self.assertRaisesRegex(SystemExit, "finite and non-negative"): + workspace.compare_fields([[1.0]], [[2.0]], + float(value) if option == "--rtol" else 0.0, + float(value) if option == "--atol" else 0.0) + + def test_invalid_initial_and_expected_timestep(self): + for value in ("nan", "inf", "-inf", "0", "-1", "1e-100", "1e100", "invalid"): + for location in ("initial", "expected"): + with self.subTest(location=location, value=value): + options = ({"initial": snapshot(mode="raw", dt=value)} if location == "initial" + else {"extra": (f"--dt={value}",)}) + self.check_driver(**options, error="timestep", launched=False) + + def test_invalid_configuration_options(self): + for extra, error in ( + (("--cells", "0,11"), "cells >= 9"), (("--cells", "13,8"), "cells >= 9"), + (("--cells", "13,11,9"), "cells >= 9"), (("--cells", "abc"), "must be nx"), + (("--steps", "0"), "steps must be positive"), (("--steps", "-1"), "steps must be positive"), + ): + self.check_driver(extra=extra, error=error, launched=False) + self.check_driver(initial=snapshot(mode="raw", steps=0), + error="positive step count", launched=False) + + def test_nonfinite_raw_fields(self): + count = 15 * 13 + for value in (math.nan, math.inf, -math.inf): + for field in ("initial", "reference", "candidate"): + with self.subTest(value=value, field=field): + self.check_driver(**{field: snapshot(mode="raw" if field == "initial" else "raw-final", + bad_value=(16, value))}, + error="non-finite", launched=field == "candidate") + for offset in (count, count + 14 * 13): + self.check_driver(initial=snapshot(mode="raw", bad_value=(offset, value)), + error="non-finite raw velocity", launched=False) + payload = snapshot(bad_value=(16, value)).split(b"\n", 3)[3] + self.check_driver(native=payload, error="non-finite", launched=False) + + def test_modes_and_payload_sizes(self): + for field in ("initial", "reference", "candidate"): + data = snapshot(mode="raw" if field == "initial" else "raw-final") + wrong_mode = data.replace(b"mode raw\n", b"mode raw-final\n") if field == "initial" \ + else data.replace(b"mode raw-final\n", b"mode raw\n") + self.check_driver(**{field: wrong_mode}, error="expected mode", launched=field == "candidate") + for payload in (data[:-1], data + b"\0" * 4): + self.check_driver(**{field: payload}, error="payload size", launched=field == "candidate") + + def test_invalid_headers(self): + for data, error in ( + (snapshot(mode="raw").replace(b"FADV2D1", b"INVALID"), "invalid header"), + (snapshot(mode="raw").replace(b"13 11 3", b"x 11 3"), "invalid dimensions"), + (snapshot(mode="raw").replace(b"13 11 3", b"13 11"), "invalid dimensions"), + ): + self.check_driver(initial=data, error=error, launched=False) + + def test_field_mismatch_boundaries_and_evolution(self): + self.check_driver(candidate=snapshot(bad_value=(16, 2.0)), error="field mismatch") + self.check_driver(candidate=snapshot(bad_value=(0, 2.0)), error="fixed boundary changed") + stationary = snapshot(mode="raw").split(b"\n", 3)[3][:15 * 13 * 4] + final = b"FADV2D1\n13 11 3 0.0001\nmode raw-final\n" + stationary + self.check_driver(reference=final, candidate=final) + self.check_driver(candidate=final, extra=("--evolution",), error="field did not change") + self.check_driver(reference=final, extra=("--evolution",), + error="field did not change", launched=False) + + def test_missing_output_cannot_reuse_stale_field(self): + (self.scratch / "workspace.final.raw").write_bytes(snapshot()) + self.check_driver(write_candidate=False, error="No such file") + + def test_shared_field_comparator_rejects_invalid_arrays(self): + for fields, error in ( + (([], []), "empty field"), (([[]], [[]]), "empty field"), + (([[1.0], [1.0, 2.0]], [[1.0], [1.0]]), "shapes differ"), + (([[math.nan]], [[1.0]]), "non-finite"), + (([[1.0]], [[math.inf]]), "non-finite"), + ): + for evolution in (False, True): + with self.assertRaisesRegex(SystemExit, error): + csv_fields.report_comparison(*fields, 1e-6, 1e-8, evolution) + for value in (math.nan, math.inf, -1.0): + with self.assertRaisesRegex(SystemExit, "finite and non-negative"): + csv_fields.report_comparison([[1.0]], [[2.0]], value, 0.0, True) + + +def check_input_guards(executable: Path): + executable = executable.resolve(strict=True) + checks = 0 + + def check(arguments, error=None): + nonlocal checks + result = subprocess.run([str(executable), *map(str, arguments)], capture_output=True, + text=True, timeout=5) + output = result.stdout + result.stderr + if error is None: + passed = result.returncode == 0 and "input validation passed" in output + else: + passed = result.returncode != 0 and error in output + if not passed: + raise SystemExit(f"guard self-test {arguments}: exit {result.returncode}\n{output}") + checks += 1 + print(f"PASS guard {arguments}: exit {result.returncode}; {error or 'accepted'}") + + def working_set(nx, ny): + return 12 * ((nx + 2) * (ny + 2) + (nx + 1) * (ny + 2) + (nx + 2) * (ny + 1)) + + cap = 512 * 1024 * 1024 + for nx, ny in ((1, 1), (5, 5), (9, 11), (51, 51), (257, 257), (2048, 2048)): + check(("working-set", nx, ny)) + square = math.isqrt(cap // 36) + while working_set(square, square) > cap: + square -= 1 + for nx, ny in ((square - 1, square - 1), (square, square), (square + 1, square + 1), (4000, 4000)): + check(("working-set", nx, ny), "exceeds 512 MiB" if working_set(nx, ny) > cap else None) + for ny in (1, 5, 9, 127, 2048): + limit = (cap // 12 - 5 * ny - 8) // (3 * ny + 5) + for nx in (limit - 1, limit, limit + 1): + error = "exceeds 512 MiB" if working_set(nx, ny) > cap else None + check(("working-set", nx, ny), error) + check(("working-set", ny, nx), error) + for nx, ny, error in ( + (2**31 - 3, 2**31 - 3, "scalar buffers size overflow"), + (800000000, 800000000, "working set size overflow"), + (2**31 - 3, 9, "exceeds 512 MiB"), (9, 2**31 - 3, "exceeds 512 MiB"), + (2**31 - 2, 9, "halo size overflow"), (9, 2**31 - 2, "halo size overflow"), + (2**31 - 1, 9, "halo size overflow"), (9, 2**31 - 1, "halo size overflow"), + (0, 9, "cells must be positive"), (9, 0, "cells must be positive"), + (-1, 9, "cells must be positive"), (9, -2**31, "cells must be positive"), + ): + check(("working-set", nx, ny), error) + for dt in ("nan", "inf", "+inf", "-inf", "0", "-0", "-1", "1e-100"): + check(("timestep", dt), "dt must be finite and positive") + for dt in ("1e-4", "1e-20", "1e-38", "1", "3.402823466e38"): + check(("timestep", dt)) + check(("defaults",)) + print(f"self-test ok: {checks} allocation-free host guard checks; no solver or GPU calls") + + +def run_self_tests(guard_executable=None): + WorkspaceValidationTests.cases = 0 + result = unittest.TextTestRunner(verbosity=2).run( + unittest.defaultTestLoader.loadTestsFromTestCase(WorkspaceValidationTests)) + if not result.wasSuccessful(): + raise SystemExit("workspace validator self-tests failed") + print(f"self-test ok: {WorkspaceValidationTests.cases} synthetic driver cases; no solver or GPU calls") + if guard_executable is not None: + check_input_guards(guard_executable) + + +if __name__ == "__main__": + run_self_tests() diff --git a/example/scripts/validate_2d_advection_workspace.py b/example/scripts/validate_2d_advection_workspace.py index a6ace19..d63bc20 100644 --- a/example/scripts/validate_2d_advection_workspace.py +++ b/example/scripts/validate_2d_advection_workspace.py @@ -2,89 +2,18 @@ from __future__ import annotations import argparse -import math -import struct import subprocess import sys from pathlib import Path -from compare_2d_fields import run_comparison as compare_csv - - -def parse_col_major_payload(payload: bytes, nx: int, ny: int): - if nx < 1 or ny < 1: - raise SystemExit("raw field dimensions must be positive") - count = (nx + 2) * (ny + 2) - if len(payload) != count * 4: - raise SystemExit(f"unexpected raw payload size: expected {count * 4} bytes") - values = struct.unpack(f"<{count}f", payload) - if not all(math.isfinite(value) for value in values): - raise SystemExit("non-finite raw scalar field") - return [[values[i + j * (nx + 2)] for j in range(ny + 2)] for i in range(nx + 2)] - - -def parse_raw_headered(path: Path, *, expected_mode=None): - with path.open("rb") as handle: - header = handle.readline().decode("ascii").strip() - if header != "FADV2D1": - raise SystemExit(f"{path}: invalid header") - dims = handle.readline().decode("ascii").strip().split() - if len(dims) != 4: - raise SystemExit(f"{path}: invalid dimensions line") - nx, ny, steps = map(int, dims[:3]) - try: - dt = struct.unpack(" atol + rtol * scale: - raise SystemExit( - f"raw field mismatch at ({i},{j}): abs diff {diff:.16e} exceeds " - f"atol+rtol*scale = {(atol + rtol * scale):.16e}" - ) - print(f"raw_max_abs_diff={max_abs:.16e}") - print(f"raw_max_rel_diff={max_rel:.16e}") - print(f"raw_changed={str(changed).lower()}") +from compare_2d_fields import ( + compare_fields, load_and_validate, report_comparison, run_comparison as compare_csv, + validate_tolerances, +) +from compare_2d_raw_fields import ( + check_final_field, check_initial_configuration, fp32_timestep, + parse_col_major_payload, parse_raw_headered, parse_raw_headerless, read_matching_final, +) def run(command, cwd: Path): @@ -93,44 +22,87 @@ def run(command, cwd: Path): def main(argv): parser = argparse.ArgumentParser(description="Validate the 2D advection-diffusion example") - parser.add_argument("--executable", required=True) - parser.add_argument("--gpu-runner", nargs="+", required=True, help="command prefix that owns the executable") + parser.add_argument("--executable") + parser.add_argument("--gpu-runner", nargs="+", help="command prefix that owns the executable") parser.add_argument("--workdir", default="build/formal-gpu/workspace-validation") - parser.add_argument("--cells", default="9,11") - parser.add_argument("--steps", default="3") - parser.add_argument("--dt", default="1e-4") + parser.add_argument("--cells", help="nx[,ny]; expected metadata with --raw-input") + parser.add_argument("--steps", type=int, help="positive step count; expected metadata with --raw-input") + parser.add_argument("--dt", help="positive FP32 timestep; expected metadata with --raw-input") parser.add_argument("--rtol", type=float, default=1e-6) parser.add_argument("--atol", type=float, default=1e-8) parser.add_argument("--raw-input") parser.add_argument("--reference-raw") parser.add_argument("--native-final-raw") parser.add_argument("--gpu-profile", action="store_true") + parser.add_argument("--evolution", action="store_true", help="require evolution beyond comparison tolerances") + parser.add_argument("--self-test", action="store_true", help="CPU-only validation helper/driver tests") + parser.add_argument("--guard-executable", type=Path, + help="also self-test the allocation-free FPM validate-advection-inputs host example") args = parser.parse_args(argv[1:]) - if args.rtol < 0.0 or args.atol < 0.0: - raise SystemExit("--rtol and --atol must be non-negative") - - workdir = Path(args.workdir) + validate_tolerances(args.rtol, args.atol) + if args.self_test: + from test_advection_validation import run_self_tests + run_self_tests(args.guard_executable) + return 0 + if args.guard_executable: + parser.error("--guard-executable requires --self-test") + if not args.executable or not args.gpu_runner: + parser.error("--executable and --gpu-runner are required unless --self-test is used") + if not args.raw_input and (args.reference_raw or args.native_final_raw): + parser.error("--reference-raw and --native-final-raw require --raw-input") + cells = None + if args.cells is not None: + try: + cells = tuple(map(int, args.cells.split(","))) + except ValueError as exc: + raise SystemExit("--cells must be nx[,ny]") from exc + if len(cells) == 1: + cells = (cells[0], cells[0]) + if len(cells) != 2 or min(cells) < 9: + raise SystemExit("--cells requires order-4 cells >= 9") + if args.steps is not None and args.steps <= 0: + raise SystemExit("--steps must be positive") + if args.dt is not None: + fp32_timestep(args.dt, "--dt") + + workdir = Path(args.workdir).resolve() workdir.mkdir(parents=True, exist_ok=True) exe = Path(args.executable).resolve() repo_root = Path(__file__).resolve().parents[2] runner = list(args.gpu_runner) if args.raw_input: + input_path = Path(args.raw_input).resolve() + initial_snapshot = parse_raw_headered(input_path, expected_mode="raw") + nx, ny, steps, _, initial = initial_snapshot + check_initial_configuration(initial_snapshot, cells if cells is not None else (nx, ny), + args.steps if args.steps is not None else steps, args.dt) + references = [] + if args.reference_raw: + references.append(("reference", read_matching_final( + Path(args.reference_raw), initial_snapshot, "reference final"))) + if args.native_final_raw: + references.append(("native", parse_raw_headerless(Path(args.native_final_raw), nx, ny))) + for label, reference in references: + check_final_field(initial, reference, label, args.rtol, args.atol, args.evolution) final_raw = workdir / "workspace.final.raw" - command = runner + [str(exe), "--solver", "workspace", "--input", args.raw_input, "--no-csv", "--dump-final-raw", str(final_raw)] + final_raw.unlink(missing_ok=True) + command = runner + [str(exe), "--solver", "workspace", "--input", str(input_path), + "--no-csv", "--dump-final-raw", str(final_raw)] if args.gpu_profile: command.append("--gpu-profile") run(command, repo_root) - candidate_nx, candidate_ny, candidate_steps, candidate_dt, candidate = parse_raw_headered(final_raw) - if args.native_final_raw: - native = parse_raw_headerless(Path(args.native_final_raw), candidate_nx, candidate_ny) - compare_fields(candidate, native, args.rtol, args.atol) - if args.reference_raw: - _, _, _, _, reference = parse_raw_headered(Path(args.reference_raw)) - compare_fields(reference, candidate, args.rtol, args.atol) + candidate = read_matching_final(final_raw, initial_snapshot, "workspace final") + check_final_field(initial, candidate, "workspace", args.rtol, args.atol, args.evolution) + for label, reference in references: + print(f"{label} versus workspace:") + report_comparison(reference, candidate, args.rtol, args.atol, False) return 0 + cells = cells if cells is not None else (9, 11) + steps = args.steps if args.steps is not None else 3 + dt = args.dt if args.dt is not None else "1e-4" operators_prefix = workdir / "operators" workspace_prefix = workdir / "workspace" operators_final = operators_prefix.with_suffix(".scalar-final.csv") @@ -144,8 +116,9 @@ def main(argv): if p.exists(): p.unlink() - operators_cmd = runner + [str(exe), "--solver", "operators", "--cells", args.cells, "--steps", args.steps, "--dt", args.dt] - workspace_cmd = runner + [str(exe), "--solver", "workspace", "--cells", args.cells, "--steps", args.steps, "--dt", args.dt] + configuration = ["--cells", ",".join(map(str, cells)), "--steps", str(steps), "--dt", dt] + operators_cmd = runner + [str(exe), "--solver", "operators", *configuration] + workspace_cmd = runner + [str(exe), "--solver", "workspace", *configuration] if args.gpu_profile: operators_cmd.append("--gpu-profile") workspace_cmd.append("--gpu-profile") @@ -157,8 +130,18 @@ def main(argv): (workdir / "workspace.scalar-final.csv").write_bytes((repo_root / "example/scripts/scalar-final.csv").read_bytes()) compare_csv(operators_final, workspace_final, rtol=args.rtol, atol=args.atol, evolution=False) compare_csv(operators_initial, workspace_initial, rtol=args.rtol, atol=args.atol, evolution=False) + if args.evolution: + for label, initial_path, final_path in ( + ("operators", operators_initial, operators_final), + ("workspace", workspace_initial, workspace_final), + ): + initial, final = load_and_validate(initial_path, final_path) + check_final_field(initial, final, label, args.rtol, args.atol, True) return 0 if __name__ == "__main__": - raise SystemExit(main(sys.argv)) + try: + raise SystemExit(main(sys.argv)) + except (OSError, ValueError) as exc: + raise SystemExit(str(exc)) from exc From ef4fb6e1cfb6eed69c588e1e92fbf739ebcb09bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Sun, 18 Oct 2026 11:20:37 -0600 Subject: [PATCH 18/18] Document input-guard checks and solve-timing limits --- benchmarks/2d-advection-metal/PERFORMANCE.md | 72 +++++++++++++++++++- benchmarks/2d-advection-metal/README.md | 14 ++++ 2 files changed, 85 insertions(+), 1 deletion(-) diff --git a/benchmarks/2d-advection-metal/PERFORMANCE.md b/benchmarks/2d-advection-metal/PERFORMANCE.md index edf3a54..7457a21 100644 --- a/benchmarks/2d-advection-metal/PERFORMANCE.md +++ b/benchmarks/2d-advection-metal/PERFORMANCE.md @@ -75,6 +75,16 @@ Check that Formal reports `gpu_timing_available=T` and positive `tg=128,1`, and positive `gpu_busy`. These solve timings exclude compilation, initialization, warmup/JIT, and field I/O. Neither command writes CSV. +**Timing caveat:** `solve_seconds` / native `solve_wall` are solve-region +wall times, not total process wall time or end-to-end benchmark time. +Formal's `final_io_seconds` measures only the analytic final-CSV path +(including scalar reconstruction when needed). The `--dump-final-raw` write +occurs after that interval; raw-input runs report `final_io_seconds=0` +even when writing a final raw field. That counter must not be used as total +final-output I/O time or added to the solve time to claim process wall time. +The raw-output timing refactor is a follow-up, not part of this input-guard +change; all recorded solve measurements retain their original scope. + For a bounded non-square command-wiring check, use `nx=17`, `ny=13`, and `steps=9` above. That is a reduced check, not a reproduction of the historical timings. @@ -120,6 +130,56 @@ cat benchmarks/2d-advection-metal/build/compare-selftest.log These synthetic comparator tests are not solver validation. +### Host input and validator checks + +The Fortran workspace estimate remains **512 MiB for three scalar and +three buffers of each face shape, all FP32**; it is not a whole-process RSS +estimate. Halo bounds, every element/byte product, and the byte-count sum +are checked before allocation, including the analytic constructor path. +Timesteps must be finite and positive; no new stability or CFL bound is +imposed. + +`example/validate-advection-inputs.f90` is an automatically discovered FPM +example used only as a host validation harness. It calls the real guards +but never allocates numerical fields, advances a timestep, or launches a +GPU kernel. Its tests cover normal grids, both sides of the working-set +boundary, halo/product/sum overflow, NaN/infinite/nonpositive timesteps, and +valid finite timesteps. **Do not run huge dimensions through the solver or +a GPU to test rejection.** + +Use the [CPU build recipe below](#controlled-cpu-measurements), replacing +`--example 2D-advection-diffusion` with `--example validate-advection-inputs` +and the copied executable with `"$cpu_build/input-guards"`. Keep the same +managed tools and CPU flags; do not use `--gpu=metal`. Then run: + +```bash +pixi run --locked --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ + "$LFORTRAN_ROOT/.pixi/envs/llvm11/bin/python" example/scripts/validate_2d_advection_workspace.py \ + --self-test --guard-executable "$cpu_build/input-guards" \ + > "$cpu_build/input-validation.log" 2>&1 +pixi run --locked --manifest-path "$LFORTRAN_ROOT/pixi.toml" -e llvm11 \ + "$LFORTRAN_ROOT/.pixi/envs/llvm11/bin/python" -O example/scripts/validate_2d_advection_workspace.py \ + --self-test --guard-executable "$cpu_build/input-guards" \ + > "$cpu_build/input-validation-optimized.log" 2>&1 +cat "$cpu_build/input-validation.log" "$cpu_build/input-validation-optimized.log" +``` + +Without `--guard-executable`, `--self-test` runs only the synthetic Python +helper/driver tests. Those deliberately mock the solver launch; they are +not evidence of CPU/GPU numerical execution. Repeat the host checks with +normal and `--fast` LFortran builds, or GFortran normal/release builds +(the measured release flags do not include `-ffast-math`). + +For real `--raw-input` validation, the initial file must have `mode raw`; +candidate/reference files must have `mode raw-final` and matching grid, +steps and FP32 timestep. Decimal spellings that round to the same FP32 +timestep are equivalent. Optional `--cells`, `--steps` and `--dt` assert +expected input metadata; otherwise the raw header is authoritative. +Nonfinite scalar/velocity fields, malformed payloads and nonfinite or +negative tolerances fail explicitly. Boundaries must remain fixed; +`--evolution` additionally requires meaningful evolution beyond the unchanged +`rtol=1e-6`, `atol=1e-8`. The validator and raw comparator share these checks. + ## Controlled CPU measurements Measured on **18 October 2026, Apple M4 Max, 64 GiB**, macOS 15.8 @@ -187,7 +247,7 @@ allocation/copy is still inside the solve. observed sample spread; these measurements establish neither a meaningful CPU benefit nor a CPU regression. The separate folded-operators/workspace ratios were **2.9067×** and **7.1055×**, respectively, on these cases. -Coefficient folding is CPU-neutral within this sample spread; the workspace +Coefficient-only differences remain unresolved within this sample spread; the workspace ratios compare workspace against **folded operators**, not old versus folded coefficients. Do not compare absolute compiler speeds or CPU/GPU times across their @@ -257,6 +317,16 @@ both managed Fortran compilers, including meaningful evolution and exact boundaries. These checks were **CPU-only**; the following section records the separate, subsequently completed GPU executions. +The subsequent host-guard follow-up passes 66 allocation-free guard checks +in each of GFortran normal/release and clean LFortran normal/`--fast`, +plus 84 synthetic driver cases in normal Python and `-O`. Actual small +Gaussian and nonzero-boundary workspace fields still pass the raw oracle +checks, and the operator/workspace CSV comparisons pass. The clean +LFortran FPM driver passes **32/32** in normal and `--fast` modes. +GFortran release passes **31/32**, with the same pre-existing 3D vector +dot-product failure reproduced using unchanged baseline library/test +sources and flags; no 3D-library or compiler workaround is included. + ## Clean-main paired GPU results Measured on **18 October 2026**, using Formal source/tooling head diff --git a/benchmarks/2d-advection-metal/README.md b/benchmarks/2d-advection-metal/README.md index 2124e0c..1715ad7 100644 --- a/benchmarks/2d-advection-metal/README.md +++ b/benchmarks/2d-advection-metal/README.md @@ -104,6 +104,20 @@ clean-main results, not merely a successful Metal link. See the [controlled CPU measurements](./PERFORMANCE.md#controlled-cpu-measurements) for the measured, workload-specific coefficient and workspace effects. +These are **solve wall times, not process wall times**. In the Formal +example, `final_io_seconds` covers the analytic final-CSV path, including +reconstruction when needed; it does **not** include `--dump-final-raw`. +Raw-input runs report that counter as zero even when writing a raw final +field. See the [timing caveat](./PERFORMANCE.md#reproduction); no end-to-end +timing claim is made. + +The [host validation checks](./PERFORMANCE.md#host-input-and-validator-checks) +exercise the Fortran working-set and finite-timestep guards without field +allocation or GPU execution. The retained workspace validation script +shares the strict raw parser, FP32 metadata, tolerance, and field checks +with the raw comparator. Huge-grid rejection tests belong only in that +allocation-free harness, never in a solver or GPU run. + ## Known Fortran compiler limitations The standalone C/Metal reference is independent of the Fortran compiler.