Scientific Computing Case Studies

Integrate with leapfrog/velocity-Verlet; extend with Barnes–Hut trees, MPI, or GPUs.

38.1 N-body gravitational simulation

Integrate with leapfrog/velocity-Verlet; extend with Barnes–Hut trees, MPI, or GPUs.

subroutine accelerations(pos, mass, acc)
  real(dp), intent(in) :: pos(:, :), mass(:)       ! pos(3, n)
  real(dp), intent(out) :: acc(:, :)
  real(dp) :: r(3), d2, inv3
  real(dp), parameter :: G = 6.674e-11_dp, soft = 1.0e-3_dp
  integer :: i, j, n
  n = size(mass)
  acc = 0.0_dp
  !$omp parallel do private(j, r, d2, inv3) schedule(dynamic)
  do i = 1, n
    do j = 1, n
      if (i == j) cycle
      r = pos(:, j) - pos(:, i)
      d2 = sum(r**2) + soft**2
      inv3 = 1.0_dp / (d2 * sqrt(d2))
      acc(:, i) = acc(:, i) + G * mass(j) * r * inv3
    end do
  end do
  !$omp end parallel do
end subroutine accelerations

38.2 Molecular dynamics (Lennard-Jones)

Neighbor lists, periodic boundary conditions (modulo/anint), velocity-Verlet, thermostats, energy conservation checks.

38.3 Computational fluid dynamics

Finite-volume/finite-difference schemes on structured grids, stencil kernels with column-major loops, halo exchange with MPI, implicit solvers via PETSc/hypre.

38.4 Climate & weather models

Large modular codebases (hundreds of thousands of lines), namelist configuration, NetCDF output, MPI + OpenMP + GPU directives, single/mixed precision for performance.

38.5 Quantum chemistry & materials

Dense linear algebra (BLAS/LAPACK/ScaLAPACK), FFTs, integral evaluation, eigen-solvers; codes like Quantum ESPRESSO, CP2K, and VASP illustrate long-lived Fortran design.

38.6 Monte Carlo simulation

program mc_pi
  use, intrinsic :: iso_fortran_env, only: dp => real64, int64
  implicit none
  integer(int64), parameter :: n = 100000000_int64
  integer(int64) :: i, hits
  real(dp) :: xy(2)
  call random_init(repeatable=.true., image_distinct=.true.)
  hits = 0
  do i = 1, n
    call random_number(xy)
    if (sum(xy**2) <= 1.0_dp) hits = hits + 1
  end do
  print '(a, f12.8)', "pi ~ ", 4.0_dp * real(hits, dp) / real(n, dp)
end program mc_pi