Case Study 2: Exposing the Heat Solver to a C Driver

"A kernel that can only be called from its own language is a program. A kernel that any language can call is a library."

Executive Summary

Where the first case study wrapped someone else's C so Fortran could use it, this one runs the traffic the other way and one step further: you will design and build a bridge so that a program written in C (and, in Chapter 15, in Python) can drive your Fortran heat solver — configure a run, execute it, and read the field back. This is the shape of real scientific software: a fast Fortran numerical core with a higher-level driver orchestrating it. Doing it well means designing three things deliberately — an interoperable configuration struct, an array layout contract that survives Fortran's column-major and C's row-major conventions without a copy, and string handling for an output filename. We build the bridge, drive it from C, and verify a hand-computed result.

Skills applied: designing an interoperable derived type to mirror a C struct (§14.5); the by-value/by-reference decision for a bridge routine (§14.3–§14.4); the column-major/row-major layout contract (§14.5); null-terminated string handling at the boundary (§14.5); and the project's bind(c) step_c shim (Project Checkpoint). It advances the heat-solver anchor: after this, the solver is callable from outside Fortran.

Background

Your solver's time step exists as the pure-Fortran step(field, alpha, dt) and, from this chapter's Project Checkpoint, the C-callable shim step_c. The new requirement: a C program should be able to set up a grid, choose the physics and the number of steps, run the simulation, and inspect the result — never touching Fortran source. You will expose one entry point, run_c, that a C main calls.

Three design decisions carry the whole bridge, and each is a §14.5 idea made concrete:

1. How does C hand Fortran the many run parameters? As one interoperable struct, not a dozen loose arguments.
2. How do the two languages agree on the 2-D grid in memory without an expensive transpose?
3. How does C pass an output filename — a null-terminated string — that Fortran can use?

Phase 1 — Design the Data Contract

Bundle the physics into a struct so the interface stays small and extensible. On the C side:

struct params_t {
    int    nsteps;
    double dx, dy, alpha, dt;
};

and the interoperable Fortran mirror — same field order, same interoperable kinds, bind(c) so the padding matches automatically:

type, bind(c) :: params_t
  integer(c_int) :: nsteps
  real(c_double) :: dx, dy, alpha, dt
end type params_t

On an LP64 platform both are 40 bytes: the int nsteps occupies bytes 0–3, four bytes of padding follow so the first double starts at byte 8, and the four doubles fill bytes 8–39. You do not compute that padding — bind(c) tells the Fortran compiler to use its companion C compiler's layout rules, and c_sizeof(p) in Fortran and sizeof(struct params_t) in C both report 40. The single discipline is to keep the two declarations in lockstep: add a field to one, add it to the other, in the same place.

Now the layout contract for the grid. Fortran will view the field as u(nx, ny), column-major, first index fastest. C stores 2-D arrays row-major, so we do not let C declare double u[ny][nx] and hope. We pass a flat block of nx*ny doubles and fix the indexing convention: Fortran's u(i+1, j+1) is C's u[i + j*nx]. Both then walk memory the same way — the fast axis i is contiguous on both sides — and no transpose is ever needed.

  Contract:  Fortran u(i+1, j+1)   <->   C  u[i + j*nx]     (i = fast/contiguous axis)

  3x3 grid, hot center at Fortran u(2,2):
     flat index of u(2,2) = (2-1) + (2-1)*3 = 4     ->  C sets u[4] = 100.0

Phase 2 — Build the Fortran Bridge

The bridge is one bind(c) routine, run_c, that loops the checkpoint's step_c over nsteps. The grid dimensions travel as by-value ints (bulletproof for sizing the array dummy), the field travels by reference (never copy a grid), and the physics travels in the struct:

! heat_bridge.f90 -- expose the heat solver to a C driver (no Fortran main).
module heat_bridge
  use, intrinsic :: iso_c_binding, only: c_int, c_double
  implicit none

  type, bind(c) :: params_t
    integer(c_int) :: nsteps
    real(c_double) :: dx, dy, alpha, dt
  end type params_t

contains

  subroutine run_c(nx, ny, u, p) bind(c, name="run_c")
    integer(c_int), value :: nx, ny               ! by value: size the grid
    real(c_double), intent(inout) :: u(nx, ny)    ! by reference: C's  double *
    type(params_t), intent(in) :: p               ! by reference: C's  const params_t *
    integer :: s
    do s = 1, p%nsteps
      call step_c(nx, ny, u, p%dx, p%dy, p%alpha, p%dt)
    end do
  end subroutine run_c

  ! The Project Checkpoint's step, as a private helper (identical math).
  subroutine step_c(nx, ny, u, dx, dy, alpha, dt)
    integer(c_int), value :: nx, ny
    real(c_double), intent(inout) :: u(nx, ny)
    real(c_double), value :: dx, dy, alpha, dt
    real(c_double) :: unew(nx, ny)
    integer :: i, j
    unew = u
    do j = 2, ny - 1
      do i = 2, nx - 1
        unew(i,j) = u(i,j) + alpha*dt * (                              &
             (u(i+1,j) - 2.0_c_double*u(i,j) + u(i-1,j)) / (dx*dx) +   &
             (u(i,j+1) - 2.0_c_double*u(i,j) + u(i,j-1)) / (dy*dy) )
      end do
    end do
    u = unew
  end subroutine step_c

end module heat_bridge

Only run_c carries bind(c), because only run_c is called from C; step_c is an internal helper and needs no C linkage. There is no program unit — C will supply main.

Phase 3 — The C Driver

C owns main. It declares the struct and the bridge prototype, fills a flat grid using the layout contract, runs, and reads the center back:

/* heat_main.c -- a C program that drives the Fortran heat solver. */
#include <stdio.h>

struct params_t {                 /* mirrors the Fortran type, bind(c) */
    int    nsteps;
    double dx, dy, alpha, dt;
};

void run_c(int nx, int ny, double *u, const struct params_t *p);

int main(void) {
    int nx = 3, ny = 3;
    double u[9] = {0};            /* column-major: u[i + j*nx] = Fortran u(i+1,j+1) */
    struct params_t p = { 2, 1.0, 1.0, 1.0, 0.1 };   /* nsteps, dx, dy, alpha, dt */

    u[1 + 1*nx] = 100.0;          /* hot center = Fortran u(2,2) */
    run_c(nx, ny, u, &p);
    printf("center after %d steps = %.2f\n", p.nsteps, u[1 + 1*nx]);
    return 0;
}

The struct is passed by address (&p), matching the by-reference type(params_t), intent(in) on the Fortran side; the flat array u is a double *, matching the by-reference u(nx,ny); and nx, ny go by value.

Phase 4 — Strings at the Boundary (Design)

The natural next feature is an output filename, and a filename is a string — the boundary's third special case. C would call an extended run_and_save_c(nx, ny, u, p, fname) with a const char *fname. Because a C string literal already carries its terminating '\0', the C caller simply passes "heat.out". The Fortran side receives it as a c_char array and must find that terminator to know where the name ends:

subroutine run_and_save_c(nx, ny, u, p, fname) bind(c, name="run_and_save_c")
  integer(c_int), value :: nx, ny
  real(c_double), intent(inout) :: u(nx, ny)
  type(params_t), intent(in) :: p
  character(kind=c_char), intent(in) :: fname(*)     ! C sees  const char *
  character(len=:), allocatable :: path
  integer :: k
  path = ""
  k = 1
  do while (fname(k) /= c_null_char)                 ! stop at the C terminator
    path = path // fname(k)
    k = k + 1
  end do
  ! ... call run_c's loop, then open(newunit=..., file=path) and write u ...
end subroutine run_and_save_c

The essential point is the do while (fname(k) /= c_null_char): a C string has no length field, so the null terminator is the length. Read past it and you scribble whatever bytes follow into your filename. (Going the other direction — Fortran handing C a filename — is the mirror image: append c_null_char, as §14.5 showed with puts.) We leave the file-writing itself to the I/O you already know; the interop content is entirely in finding the terminator.

Phase 5 — Build, Run, and Verify

Compile each side, link with gfortran (Fortran is present, so libgfortran must come along):

$ gfortran -std=f2018 -Wall -c heat_bridge.f90
$ gcc -c heat_main.c
$ gfortran heat_bridge.o heat_main.o -o heat
$ ./heat
center after 2 steps = 36.00

Verify the number by hand. The grid is $3 \times 3$ with a single interior point $u(2,2)$, all else zero, $\alpha = \Delta x = \Delta y = 1$, $\Delta t = 0.1$.

• Step 1. Laplacian at the center $= (0 - 2\cdot100 + 0) + (0 - 2\cdot100 + 0) = -400$, so $u(2,2) \to 100 + 0.1(-400) = 60$.
• Step 2. Now $u(2,2) = 60$: Laplacian $= (0 - 120 + 0) + (0 - 120 + 0) = -240$, so $u(2,2) \to 60 + 0.1(-240) = 36$.

C reads the center at flat index u[1 + 1*3] = u[4], our contract's mapping of Fortran u(2,2), and prints 36.00. The two languages agreed on the struct, on the grid layout, and on the entry point — and the result is exactly what pure Fortran would have produced.

A cheap but revealing sanity check: change the C driver to set u[1 + 0*nx] = 100.0 instead (Fortran u(2,1), a boundary cell) and rerun — the center stays 0.00, confirming that the layout contract really does place your hot spot where you think it does. Break the contract instead (index the hot spot row-major, u[1*nx + 1]) and the physics quietly moves; that silent relocation is the whole reason the contract must be explicit.

Discussion Questions

1. We passed nx and ny as by-value ints and could have put them in params_t. Give one reason to keep the grid dimensions as separate scalar arguments rather than struct fields, thinking about how the array dummy u(nx, ny) is declared.
2. The layout contract keeps Fortran column-major and asks C to index u[i + j*nx]. Suppose instead you insisted C keep its natural u[j][i] row-major grid. What would the bridge have to do on every call, and what would it cost for a $4000 \times 4000$ grid (recall Exercise 14.20)?
3. run_c takes params_t by reference (const params_t *) rather than by value. Why is by-reference the right default for a struct argument, and when might by value be acceptable?

Your Turn: Extensions

• Option A. Add a double *max_change output argument to run_c so the C driver learns the largest absolute change over the whole run (a convergence signal). What interoperable declaration does it get on the Fortran side, and why can it not be value?
• Option B. Implement Phase 4 fully: write run_and_save_c, have it write the final field to the named file as text, and drive it from C with run_and_save_c(nx, ny, u, &p, "heat.out"). Confirm the file appears and its contents match the hand-computed field.
• Option C. Replace the C driver with a Python one using ctypes, calling the very same run_c symbol in a shared library built from heat_bridge.f90. You will meet the same struct and the same layout contract — a direct preview of Chapter 15, where f2py automates most of this.

Key Takeaways

• A reusable numerical kernel needs a designed boundary: an interoperable struct for configuration, an explicit array-layout contract, and a plan for strings. Each is a §14.5 idea, and together they turn a program into a callable library.
• The layout contract — Fortran u(i+1,j+1) = C u[i + j*nx] — lets a 2-D grid cross the boundary with no copy. Insisting on each language's "natural" 2-D indexing instead forces a per-call transpose you cannot afford.
• bind(c) on the one entry point C calls is enough; internal helpers stay ordinary Fortran. Keep the struct definitions on the two sides in lockstep, and let gfortran drive the link.
• Verify a mixed-language result against hand computation, and probe the layout contract with a boundary cell — because the failure mode here is not a crash but a silently relocated physics.