HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
1D Inviscid Burgers Equation - Sine Wave

Location: hypar/Examples/1D/Burgers/SineWave (This directory contains all the input files needed to run this case.)

Governing equations: 1D Burgers Equation (burgers.h)

References:

  • Ghosh, D., Baeder, J. D., "Compact Reconstruction Schemes with Weighted ENO Limiting for Hyperbolic Conservation Laws", SIAM Journal on Scientific Computing, 34 (3), 2012, A1678–A1706

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

Initial solution: \(u\left(x,0\right) = \frac{1}{2 \pi t_s}\sin\left(2\pi x\right)\), \(t_s=2\) (time to shock formation)

Exact solution: \(u\left(x,y,t\right) = \frac{1}{2 \pi t_s} \sin\left(2\pi \left(x-u t\right)\right)\) (needs to be computed iteratively)
Numerical Method:

Input files required:

solver.inp

begin
ndims 1
nvars 1
size 80
iproc 1
ghost 3
n_iter 20
time_scheme rk
time_scheme_type ssprk3
hyp_space_scheme crweno5
conservation_check yes
dt 0.1
screen_op_iter 1
file_op_iter 3
ip_file_type ascii
op_file_format text
op_overwrite no
model burgers
end

boundary.inp

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

physics.inp

begin
end

lusolver.inp (optional)

begin
reducedsolvetype jacobi
evaluate_norm 1
maxiter 10
atol 1e-12
rtol 1e-10
verbose 0
end

weno.inp (optional)

begin
mapped 1
borges 0
yc 0
no_limiting 0
epsilon 0.000001
p 2.0
rc 0.3
xi 0.001
end

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

/*
Code to generate the initial and exact
solutions for:
Case: Sine Wave
Model: Burger1D
Needs: solver.inp
Writes out: initial.inp
If the final time is less than shock
formation time, then the exact solution
is also written out.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
int main()
{
/* shock formation time */
double ts = 2.0;
/* maximum number of iterations for computing exact solution */
int MAX_ITER = 10000;
/* Convergence tolerance for computing exact solution */
double tolerance = 1e-15;
/* value of pi */
double pi = 4.0*atan(1.0);
int NI, ndims, niter;
double dt;
char ip_file_type[50];
strcpy(ip_file_type,"ascii");
FILE *in;
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, "ip_file_type")) fscanf(in,"%s",ip_file_type);
else if (!strcmp(word, "n_iter")) fscanf(in,"%d" , &niter);
else if (!strcmp(word, "dt")) fscanf(in,"%lf", &dt);
}
} 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);
int i;
double dx = 1.0 / ((double)NI);
double tf = ((double)niter) * dt;
printf("Final Time: %lf\n",tf);
/* initial solution */
{
double *x, *u;
x = (double*) calloc (NI, sizeof(double));
u = (double*) calloc (NI, sizeof(double));
for (i = 0; i < NI; i++){
x[i] = i*dx;
u[i] = sin(2.0*pi*x[i]) / (ts*2.0*pi);
}
FILE *out;
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);
}
free(x);
free(u);
}
/* exact solution */
if (tf < ts) {
double *x, *u;
x = (double*) calloc (NI, sizeof(double));
u = (double*) calloc (NI, sizeof(double));
for (i = 0; i < NI; i++){
x[i] = i*dx;
u[i] = sin(2.0*pi*x[i]) / (ts*2.0*pi);
}
int k;
printf("Computing exact solution iteratively...\n");
for (k = 0; k < MAX_ITER; k++) {
double maxres = 0;
for (i = 0; i < NI; i++){
double new_u = sin(2.0*pi*(x[i]-u[i]*tf)) / (ts*2.0*pi);
double res = sqrt((new_u-u[i])*(new_u-u[i]));
if (res > maxres) maxres = res;
u[i] = new_u;
}
printf(" iter=%6d, max res=%1.6e\n", k, maxres);
if (maxres < tolerance) break;
}
FILE *out;
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);
} else {
printf("Final time (%f) greater than shock formation time (%f).\n", tf, ts);
printf("Exact solution not available.\n");
}
return(0);
}

Note: The exact solution is available only if the final time is less than \(t_s\) above.

Output:

After running the code, there should be 8 output files op_00000.dat, op_00001.dat, ... op_00007.dat; the first one is the solution at \(t=0\) and the final one is the solution at \(t=2\). 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: The following figure is obtained by plotting op_00000.dat (initial), ..., op_00007.dat (t=2).

Solution_1DBurgersSine.png

Since HyPar::ConservationCheck is set to yes in solver.inp, the code checks for conservation error and prints it to screen, as well as the file conservation.dat:

80 1 1.0000000000000001E-01 3.2797115717686509E-18

The numbers are: number of grid points (HyPar::dim_global), number of processors (MPIVariables::iproc), time step size (HyPar::dt), and conservation error (HyPar::ConservationError).

Expected screen output:

HyPar - Parallel (MPI) version with 1 processes
Compiled with PETSc time integration.
Allocated simulation object(s).
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. : 20
Restart iteration : 0
Time integration scheme : rk (ssprk3)
Spatial discretization scheme (hyperbolic) : crweno5
Split hyperbolic flux term? : no
Interpolation type for hyperbolic term : characteristic
Spatial discretization type (parabolic ) : nonconservative-1stage
Spatial discretization scheme (parabolic ) : 2
Time Step : 1.000000E-01
Check for conservation : yes
Screen output iterations : 1
File output iterations : 3
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 : burgers
Partitioning domain and allocating data arrays.
Reading array from ASCII file initial.inp (Serial mode).
Volume integral of the initial solution:
0: -4.3368086899420177E-18
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.
Reading WENO parameters from weno.inp.
Initializing physics. Model = "burgers"
Reading physical model inputs from file "physics.inp".
Setting up time integration.
Solving in time (from 0 to 20 iterations)
Writing solution file op_00000.dat.
Iteration: 1 Time: 1.000E-01 Max CFL: 6.366E-01 Max Diff. No.: -1.000E+00 Norm: 1.4067E-03 Conservation loss: 9.7578E-19
Iteration: 2 Time: 2.000E-01 Max CFL: 6.364E-01 Max Diff. No.: -1.000E+00 Norm: 1.4084E-03 Conservation loss: 2.4395E-19
Iteration: 3 Time: 3.000E-01 Max CFL: 6.365E-01 Max Diff. No.: -1.000E+00 Norm: 1.4120E-03 Conservation loss: 2.2497E-18
Writing solution file op_00001.dat.
Iteration: 4 Time: 4.000E-01 Max CFL: 6.366E-01 Max Diff. No.: -1.000E+00 Norm: 1.4174E-03 Conservation loss: 7.1828E-19
Iteration: 5 Time: 5.000E-01 Max CFL: 6.362E-01 Max Diff. No.: -1.000E+00 Norm: 1.4247E-03 Conservation loss: 2.8867E-18
Iteration: 6 Time: 6.000E-01 Max CFL: 6.365E-01 Max Diff. No.: -1.000E+00 Norm: 1.4341E-03 Conservation loss: 6.2613E-18
Writing solution file op_00002.dat.
Iteration: 7 Time: 7.000E-01 Max CFL: 6.365E-01 Max Diff. No.: -1.000E+00 Norm: 1.4457E-03 Conservation loss: 2.9883E-18
Iteration: 8 Time: 8.000E-01 Max CFL: 6.362E-01 Max Diff. No.: -1.000E+00 Norm: 1.4597E-03 Conservation loss: 4.1335E-19
Iteration: 9 Time: 9.000E-01 Max CFL: 6.366E-01 Max Diff. No.: -1.000E+00 Norm: 1.4764E-03 Conservation loss: 1.3553E-19
Writing solution file op_00003.dat.
Iteration: 10 Time: 1.000E+00 Max CFL: 6.365E-01 Max Diff. No.: -1.000E+00 Norm: 1.4961E-03 Conservation loss: 1.3431E-17
Iteration: 11 Time: 1.100E+00 Max CFL: 6.363E-01 Max Diff. No.: -1.000E+00 Norm: 1.5195E-03 Conservation loss: 6.1664E-19
Iteration: 12 Time: 1.200E+00 Max CFL: 6.366E-01 Max Diff. No.: -1.000E+00 Norm: 1.5470E-03 Conservation loss: 4.2284E-18
Writing solution file op_00004.dat.
Iteration: 13 Time: 1.300E+00 Max CFL: 6.363E-01 Max Diff. No.: -1.000E+00 Norm: 1.5795E-03 Conservation loss: 4.8654E-18
Iteration: 14 Time: 1.400E+00 Max CFL: 6.365E-01 Max Diff. No.: -1.000E+00 Norm: 1.6183E-03 Conservation loss: 5.1500E-18
Iteration: 15 Time: 1.500E+00 Max CFL: 6.366E-01 Max Diff. No.: -1.000E+00 Norm: 1.6649E-03 Conservation loss: 2.8867E-18
Writing solution file op_00005.dat.
Iteration: 16 Time: 1.600E+00 Max CFL: 6.362E-01 Max Diff. No.: -1.000E+00 Norm: 1.7213E-03 Conservation loss: 5.4142E-18
Iteration: 17 Time: 1.700E+00 Max CFL: 6.365E-01 Max Diff. No.: -1.000E+00 Norm: 1.7870E-03 Conservation loss: 1.1520E-18
Iteration: 18 Time: 1.800E+00 Max CFL: 6.365E-01 Max Diff. No.: -1.000E+00 Norm: 1.8481E-03 Conservation loss: 1.0300E-18
Writing solution file op_00006.dat.
Iteration: 19 Time: 1.900E+00 Max CFL: 6.362E-01 Max Diff. No.: -1.000E+00 Norm: 1.8811E-03 Conservation loss: 8.0773E-18
Iteration: 20 Time: 2.000E+00 Max CFL: 6.366E-01 Max Diff. No.: -1.000E+00 Norm: 1.8334E-03 Conservation loss: 3.2797E-18
Writing solution file op_00007.dat.
Completed time integration (Final time: 2.000000).
Conservation Errors:
3.2797115717686509E-18
Solver runtime (in seconds): 2.4938999999999999E-02
Total runtime (in seconds): 3.4243000000000003E-02
Deallocating arrays.
Finished.