Laplace solver

laplace.c

// Copyright [2023] Taylor Howell
// Laplace solver

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

// set boundary conditions for grid
void set_boundaries(double* grid, int width, int height, double top,
                    double bottom, double left, double right) {
  // top + bottom
  for (int i = 0; i < width; i++) {
    grid[i] = top;
    grid[(height - 1) * width + i] = bottom;
  }

  // left + right
  for (int i = 0; i < height; i++) {
    grid[i * width] = left;
    grid[i * width + (width - 1)] = right;
  }
}

// iteratively solve grid until convergence to tolerance
void solve(double* grid, int width, int height, int max_iter,
           double tolerance) {
  // update iterations
  for (int k = 0; k < max_iter; k++) {
    // track max error for early termination
    double max_error = 0.0;

    // loop over interior of grid
    for (int i = 1; i < height - 1; i++) {
      for (int j = 1; j < width - 1; j++) {
        // average neighbor grid cells
        double ut = grid[(i - 1) * width + j];
        double ub = grid[(i + 1) * width + j];
        double ul = grid[i * width + j - 1];
        double ur = grid[i * width + j + 1];
        double gij = 0.25 * (ut + ub + ul + ur);

        // compute error
        double error = fabs(gij - grid[i * width + j]);
        if (error > max_error) {
          max_error = error;
        }

        // update grid value
        grid[i * width + j] = gij;
      }
    }

    // check convergence
    if (max_error < tolerance) {
      printf("[success after %i iterations]\n\n", k);
      return;
    }
  }

  // failed to converge
  printf("[failed after %i iterations]\n\n", max_iter);
}

// print 2D grid to terminal
void print_grid(const double* grid, int width, int height) {
  // loop over grid
  for (int i = 0; i < height; i++) {
    for (int j = 0; j < width; j++) {
      printf("%.2f ", grid[i * width + j]);
    }
    printf("\n");
  }
  printf("\n");
}

int main() {
  printf("Laplace solver\n\n");

  // dimensions
  int width = 10;
  int height = 10;

  // grid
  double* grid = malloc(width * height * sizeof(double));

  // initialize
  for (int i = 0; i < width * height; i++) {
    grid[i] = 0.0;
  }

  // initial grid
  printf("initial grid\n");
  print_grid(grid, width, height);

  // set boundary conditions
  double top = 0.0;
  double bottom = 0.0;
  double left = 0.0;
  double right = 1.0;
  set_boundaries(grid, width, height, top, bottom, left, right);

  // grid w/ boundary conditions
  printf("grid w/ boundary conditions\n");
  print_grid(grid, width, height);

  // solve
  int max_iter = 1000;
  double tolerance = 1.0e-6;
  printf("solve: ");
  solve(grid, width, height, max_iter, tolerance);

  // grid at steady state
  printf("grid at steady state\n");
  print_grid(grid, width, height);

  // free memory
  free(grid);

  return 0;
}