Location: hypar/Examples/1D/LinearDiffusion/SineWave (This directory contains all the input files needed to run this case. If there is a Run.m, run it in MATLAB to quickly set up, run, and visualize the example).

Governing equations: 1D Linear Diffusion Equation (linearadr.h)

Domain: \(0 \le x < 1\), "periodic" (_PERIODIC_) boundary conditions

Initial solution: \(u\left(x,0\right) = \sin\left(2\pi x\right)\)

Numerical Method:

Spatial discretization (parabolic): 2nd order (Interp2PrimSecondOrder()), conservative (ParabolicFunctionCons1Stage())
Time integration: SSPRK3 (TimeRK(), _RK_SSP3_)

Input files required:

solver.inp

begin
  ndims             1
  nvars             1
  size              80
  iproc             1
  ghost             3
  n_iter            200
  time_scheme       rk
  time_scheme_type  ssprk3
  par_space_type    conservative-1stage
  par_space_scheme  2
  dt                0.05
  screen_op_iter    10
  file_op_iter      20
  op_file_format    text
  op_overwrite      no
  model             linear-advection-diffusion-reaction
end

boundary.inp

2
periodic      0     1     0     0
periodic      0    -1     0     0

physics.inp

begin
  diffusion 0.001
end

To generate initial.inp (initial solution) and exact.inp (exact solution), compile and run the following code in the run directory.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
int main(){
  
  double pi = 4.0*atan(1.0);
    int NI,ndims,niter;
  double nu, dt,final_time;
  char ip_file_type[50];
  strcpy(ip_file_type,"ascii");
  FILE *in, *out;
  printf("Reading file \"solver.inp\"...\n");
  in = fopen("solver.inp","r");
  if (!in) printf("Error: Input file \"solver.inp\" not found. Default values will be used.\n");
  else {
    char word[500];
    fscanf(in,"%s",word);
    if (!strcmp(word, "begin")){
      while (strcmp(word, "end")){
        fscanf(in,"%s",word);
        if (!strcmp(word, "ndims"))       fscanf(in,"%d",&ndims);
        else if (!strcmp(word, "size"))   fscanf(in,"%d",&NI);
        else if (!strcmp(word, "dt"))     fscanf(in,"%lf",&dt);
        else if (!strcmp(word, "n_iter")) fscanf(in,"%d",&niter);
        else if (!strcmp(word, "ip_file_type")) fscanf(in,"%s",ip_file_type);
      }
    } else printf("Error: Illegal format in solver.inp. Crash and burn!\n");
  }
  fclose(in);
  if (ndims != 1) {
    printf("ndims is not 1 in solver.inp. this code is to generate 1D initial conditions\n");
    return(0);
  }
    printf("Grid:\t\t\t%d\n",NI);
  final_time = (double)niter * dt;
    printf("Final Time:\t\t\t%lf\n",final_time);
  printf("Reading file \"physics.inp\"...\n");
  in = fopen("physics.inp","r");
  if (!in) {
    printf("Error: Input file \"physics.inp\" not found. Default values will be used.\n");
    nu = 1.0;
  } else {
    char word[500];
    fscanf(in,"%s",word);
    if (!strcmp(word, "begin")){
      while (strcmp(word, "end")){
        fscanf(in,"%s",word);
        if (!strcmp(word,"diffusion")) fscanf(in,"%lf",&nu);
      }
    } else printf("Error: Illegal format in solver.inp. Crash and burn!\n");
  }
  fclose(in);
    printf("Diffusion Coeff:\t\t\t%lf\n",nu);
    int i;
    double dx = 1.0 / ((double)NI);
    double *x, *u;
    x = (double*) calloc (NI, sizeof(double));
    u = (double*) calloc (NI, sizeof(double));
  /* set grid */
    for (i = 0; i < NI; i++) x[i] = i*dx;
  /* Initial  solution */
    for (i = 0; i < NI; i++) u[i] = sin(2*pi*x[i]);
  if (!strcmp(ip_file_type,"ascii")) {
    printf("Writing ASCII initial solution file initial.inp\n");
      out = fopen("initial.inp","w");
    for (i = 0; i < NI; i++)  fprintf(out,"%lf ",x[i]);
    fprintf(out,"\n");
      for (i = 0; i < NI; i++)  fprintf(out,"%lf ",u[i]);                       
    fprintf(out,"\n");
      fclose(out);
  } else if ((!strcmp(ip_file_type,"binary")) || (!strcmp(ip_file_type,"bin"))) {
    printf("Writing binary initial solution file initial.inp\n");
    out = fopen("initial.inp","wb");
    fwrite(x,sizeof(double),NI,out);
    fwrite(u,sizeof(double),NI,out);
    fclose(out);
  }
  /* Exact solution */
    for (i = 0; i < NI; i++) u[i] = exp(-nu*4*pi*pi*final_time)*sin(2*pi*x[i]);
  if (!strcmp(ip_file_type,"ascii")) {
    printf("Writing ASCII exact solution file exact.inp\n");
      out = fopen("exact.inp","w");
    for (i = 0; i < NI; i++)  fprintf(out,"%lf ",x[i]);
    fprintf(out,"\n");
      for (i = 0; i < NI; i++)  fprintf(out,"%lf ",u[i]);                       
    fprintf(out,"\n");
      fclose(out);
  } else if ((!strcmp(ip_file_type,"binary")) || (!strcmp(ip_file_type,"bin"))) {
    printf("Writing binary exact solution file exact.inp\n");
    out = fopen("exact.inp","wb");
    fwrite(x,sizeof(double),NI,out);
    fwrite(u,sizeof(double),NI,out);
    fclose(out);
  }
    free(x);
    free(u);
    return(0);
}

Output:

After running the code, there should be 11 output files op_00000.dat, op_00001.dat, ... op_00010.dat; the first one is the solution at \(t=0\) and the final one is the solution at \(t=10\). Since HyPar::op_overwrite is set to no in solver.inp, separate files are written for solutions at each output time. All the files are ASCII text (HyPar::op_file_format is set to text in solver.inp). In these files, the first column is grid index, the second column is x-coordinate, and the third column is the solution.

Solutions at t=0,2,4,6,8,10: The following figure is obtained by plotting op_00000.dat (t=0, initial), op_00002.dat (t=2), op_00004.dat (t=4), op_00006.dat (t=6), op_00008.dat (t=8) and op_00010.dat (t=10, final).

Since the exact solution is available at the final time , the numerical errors are calculated and reported on screen (see below) as well as errors.dat:

80 1 5.0000000000000003E-02 2.0258640515128642E-04 2.0266155730691682E-04 2.0361464428477190E-04 1.2026999999999999E-02 1.2685000000000000E-02

The numbers are: number of grid points (HyPar::dim_global), number of processors (MPIVariables::iproc), time step size (HyPar::dt), L1, L2, and L-infinity errors (HyPar::error), solver wall time (seconds) (i.e., not accounting for initialization, and cleaning up), and total wall time.

Expected screen output:

HyPar - Parallel (MPI) version with 1 processes
Compiled with PETSc time integration.
Reading solver inputs from file "solver.inp".
    No. of dimensions                          : 1
    No. of variables                           : 1
    Domain size                                : 80 
    Processes along each dimension             : 1 
    No. of ghosts pts                          : 3
    No. of iter.                               : 200
    Restart iteration                          : 0
    Time integration scheme                    : rk (ssprk3)
    Spatial discretization scheme (hyperbolic) : 1
    Split hyperbolic flux term?                : no
    Interpolation type for hyperbolic term     : characteristic
    Spatial discretization type   (parabolic ) : conservative-1stage
    Spatial discretization scheme (parabolic ) : 2
    Time Step                                  : 5.000000E-02
    Check for conservation                     : no
    Screen output iterations                   : 10
    File output iterations                     : 20
    Initial solution file type                 : ascii
    Initial solution read mode                 : serial
    Solution file write mode                   : serial
    Solution file format                       : text
    Overwrite solution file                    : no
    Physical model                             : linear-advection-diffusion-reaction
Partitioning domain.
Allocating data arrays.
Reading array from ASCII file initial.inp (Serial mode).
Volume integral of the initial solution:
 0:  4.9005938196344800E-17
Reading boundary conditions from "boundary.inp".
  Boundary                       periodic:  Along dimension  0 and face +1
  Boundary                       periodic:  Along dimension  0 and face -1
2 boundary condition(s) read.
Initializing solvers.
Initializing physics. Model = "linear-advection-diffusion-reaction"
Reading physical model inputs from file "physics.inp".
Setting up time integration.
Solving in time (from 0 to 200 iterations)
Writing solution file op_00000.dat.
Iteration:      10  Time: 5.000E-01  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.3692E-03  
Iteration:      20  Time: 1.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.3424E-03  
Writing solution file op_00001.dat.
Iteration:      30  Time: 1.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.3162E-03  
Iteration:      40  Time: 2.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.2905E-03  
Writing solution file op_00002.dat.
Iteration:      50  Time: 2.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.2653E-03  
Iteration:      60  Time: 3.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.2405E-03  
Writing solution file op_00003.dat.
Iteration:      70  Time: 3.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.2163E-03  
Iteration:      80  Time: 4.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.1925E-03  
Writing solution file op_00004.dat.
Iteration:      90  Time: 4.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.1692E-03  
Iteration:     100  Time: 5.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.1464E-03  
Writing solution file op_00005.dat.
Iteration:     110  Time: 5.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.1240E-03  
Iteration:     120  Time: 6.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.1020E-03  
Writing solution file op_00006.dat.
Iteration:     130  Time: 6.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.0805E-03  
Iteration:     140  Time: 7.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.0594E-03  
Writing solution file op_00007.dat.
Iteration:     150  Time: 7.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.0387E-03  
Iteration:     160  Time: 8.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 1.0184E-03  
Writing solution file op_00008.dat.
Iteration:     170  Time: 8.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 9.9853E-04  
Iteration:     180  Time: 9.000E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 9.7902E-04  
Writing solution file op_00009.dat.
Iteration:     190  Time: 9.500E+00  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 9.5989E-04  
Iteration:     200  Time: 1.000E+01  Max CFL: 0.000E+00  Max Diff. No.: 3.200E-01  Norm: 9.4114E-04  
Writing solution file op_00010.dat.
Completed time integration (Final time: 10.000000).
Reading array from ASCII file exact.inp (Serial mode).
Computed errors:
  L1         Error           : 2.0258640515128642E-04
  L2         Error           : 2.0266155730691682E-04
  Linfinity  Error           : 2.0361464428477190E-04
Conservation Errors:
    0.0000000000000000E+00
Solver runtime (in seconds): 1.2026999999999999E-02
Total  runtime (in seconds): 1.2685000000000000E-02
Deallocating arrays.
Finished.