Solve the compressible Navier-Stokes equation in 1d with a polytropic equation of state.

isentropic.c

/*
 * AUTHOR: Jonathan Zrake
 *
 * DATE: 3/31/2013
 * 
 * PURPOSE: Solve the compressible Navier-Stokes equation in 1d with a
 * polytropic equation of state.
 *
 * SCHEME: semi-implicit Runge-Kutta / Backward Euler
 *
 *   The hyperbolic terms (Euler equations) are stepped explicitly using a
 *   4th-order centered finite difference stencil.
 *
 *   The parabolic term controlling viscous dissipation is evolved implicitly in
 *   order not to dominate the CFL constraint.
 *
 *   Performance is always stable for CFL < 1 or even a little larger. If the
 *   viscosity is too low and a shock is permitted to form, then the finite
 *   differences will produce over-shoots, but not crashes.
 *
 */

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

#define NX 1024
#define PI (4*atan(1))

static double CFL = 0.9;
static double Gamma = 1.4;
static double Kappa = 1.0;
static double Viscosity = 1e-3;
static void diffusion(double *u, double dt);

double x_at(int i)
{
  return (double)i / NX - 0.5;
}
double exact_solution_u(double x, double t)
{
  return 0.0;
}
double exact_solution_d(double x, double t)
{
  return 1.0 + 0.2 * exp(-x*x/0.002);
}
double diff4(double *u, int i)
{
  int ip1 = i + 1; if (ip1 >= NX) ip1 -= NX;
  int im1 = i - 1; if (im1 <   0) im1 += NX;
  int ip2 = i + 2; if (ip2 >= NX) ip2 -= NX;
  int im2 = i - 2; if (im2 <   0) im2 += NX;
  return u[im2]/12 - 2*u[im1]/3 + 2*u[ip1]/3 - u[ip2]/12;
}

int main()
{
  double x;
  double u0[NX], u1[NX], u2[NX];
  double d0[NX], d1[NX], d2[NX];
  double dt, dx = 1.0 / NX;
  double t = 0.0;
  double cs2, max_cs2;
  double dudx, dudt;
  double dddx, dddt;

  int iter = 0;
  int i;

  for (i=0; i<NX; ++i) {
    x = x_at(i);
    u0[i] = exact_solution_u(x, t);
    d0[i] = exact_solution_d(x, t);
  }

  while (iter < 8*NX) {

    if (iter % 16 == 0) {
      char fname[128];
      snprintf(fname, 128, "data/out-%05d.dat", iter);
      FILE *outf = fopen(fname, "w");

      for (i=0; i<NX; ++i) {
        x = x_at(i);
        fprintf(outf, "%+12.10e %+12.10e %+12.10e %+12.10e %+12.10e\n",
                x_at(i), d0[i], u0[i],
                exact_solution_d(x, t),
                exact_solution_u(x, t));
      }
      fclose(outf);
    }

    max_cs2 = 0.0;
    for (i=0; i<NX; ++i) {
      cs2 = Kappa * Gamma * pow(d0[i], Gamma - 1);
      if (max_cs2 < cs2) max_cs2 = cs2;
    }
    dt = CFL * dx / sqrt(max_cs2);

    diffusion(u0, dt);
    for (i=0; i<NX; ++i) {
      cs2  = Kappa * Gamma * pow(d0[i], Gamma - 1);
      dudx = diff4(u0, i) / dx;
      dddx = diff4(d0, i) / dx;
      dddt =       -u0[i] * dddx - d0[i] * dudx;
      dudt = -cs2 / d0[i] * dddx - u0[i] * dudx;
      u2[i] = u0[i] + dt * dudt;
      d2[i] = d0[i] + dt * dddt;
    }
    for (i=0; i<NX; ++i) u1[i] = u2[i];
    for (i=0; i<NX; ++i) d1[i] = d2[i];

    diffusion(u1, dt);
    for (i=0; i<NX; ++i) {
      cs2  = Kappa * Gamma * pow(d0[i], Gamma - 1);
      dudx = diff4(u1, i) / dx;
      dddx = diff4(d1, i) / dx;
      dddt =       -u1[i] * dddx - d1[i] * dudx;
      dudt = -cs2 / d1[i] * dddx - u1[i] * dudx;
      u2[i] = (3./4.) * u0[i] + (1./4.) * (u1[i] + dt * dudt);
      d2[i] = (3./4.) * d0[i] + (1./4.) * (d1[i] + dt * dddt);
    }
    for (i=0; i<NX; ++i) u1[i] = u2[i];
    for (i=0; i<NX; ++i) d1[i] = d2[i];

    diffusion(u1, dt);
    for (i=0; i<NX; ++i) {
      cs2  = Kappa * Gamma * pow(d0[i], Gamma - 1);
      dudx = diff4(u1, i) / dx;
      dddx = diff4(d1, i) / dx;
      dddt =       -u1[i] * dddx - d1[i] * dudx;
      dudt = -cs2 / d1[i] * dddx - u1[i] * dudx;
      u2[i] = (1./3.) * u0[i] + (2./3.) * (u1[i] + dt * dudt);
      d2[i] = (1./3.) * d0[i] + (2./3.) * (d1[i] + dt * dddt);
    }
    for (i=0; i<NX; ++i) u0[i] = u2[i];
    for (i=0; i<NX; ++i) d0[i] = d2[i];

    iter += 1;
    t += dt;
    printf("t=%f dt=%f\n", t, dt);
  }

  return 0;
}


void diffusion(double *u, double dt)
/*
 *
 * Take a diffusive Backward Euler step, using the DFT to work in the Fourier
 * domain.
 *
 *              u^{n+1} = u^{n} + dt * nu \nabla^2 u^{n+1}
 *
 */
{
  int i;
  double k;
  fftw_complex *ux = (fftw_complex*) fftw_malloc(NX * sizeof(fftw_complex));
  fftw_complex *uk = (fftw_complex*) fftw_malloc(NX * sizeof(fftw_complex));

  fftw_plan fwd = fftw_plan_dft_1d(NX, ux, uk, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_plan rev = fftw_plan_dft_1d(NX, uk, ux, FFTW_BACKWARD, FFTW_ESTIMATE);

  for (i=0; i<NX; ++i) {
    ux[i][0] = u[i];
    ux[i][1] = 0.0;
  }
  fftw_execute(fwd);

  /* For compressible NS equations in 1d, the effective viscosity is larger by a
     factor of 4/3. */
  for (i=0; i<NX; ++i) {
    k = 2 * PI * (i < NX/2 ? i : i-NX);
    uk[i][0] /= (1 + dt * 4./3 * Viscosity * k*k);
    uk[i][1] /= (1 + dt * 4./3 * Viscosity * k*k);
  }
  fftw_execute(rev);

  for (i=0; i<NX; ++i) {
    u[i] = ux[i][0] / NX; // divide by NX to correct for normalization
  }

  fftw_free(ux);
  fftw_free(uk);
  fftw_destroy_plan(fwd);
  fftw_destroy_plan(rev);
}