Location: hypar/Examples/2D/Burgers/SineWave (This directory contains all the input files needed to run this case.)
Governing equations: 2D Burgers Equation (burgers.h)
Domain: \(1 \le x,y < 1\), "periodic" (_PERIODIC_) boundary conditions on all boundaries.
Initial solution: \(u\left(x,y,0\right) = \frac{1}{2 \pi t_s} \sin\left(2\pi x\right) \sin\left(2\pi y\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) \sin\left(2\pi \left(y-u t\right)\right)\) (needs to be computed iteratively)
Numerical Method:
Input files required:
solver.inp
begin
ndims 2
nvars 1
size 64 64
iproc 2 2
ghost 3
n_iter 16
restart_iter 0
time_scheme rk
time_scheme_type 44
hyp_space_scheme crweno5
dt 0.1
conservation_check yes
screen_op_iter 1
file_op_iter 4
op_file_format tecplot2d
ip_file_type binary
input_mode serial
output_mode serial
op_overwrite no
model burgers
end
boundary.inp
4
periodic 0 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+00
periodic 0 -1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+00
periodic 1 1 0.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
periodic 1 -1 0.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
physics.inp (specifies \(a_x\) and \(a_y\))
lusolver.inp (optional)
begin
reducedsolvetype jacobi
evaluate_norm 1
maxiter 10
atol 9.9999999999999998e-13
rtol 1.0000000000000000e-10
verbose 0
end
weno.inp (optional)
begin
mapped 0
borges 0
yc 1
no_limiting 0
epsilon 9.9999999999999995e-07
p 2.0000000000000000e+00
rc 2.9999999999999999e-01
xi 1.0000000000000000e-03
tol 9.9999999999999998e-17
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>
{
double ts = 2.0;
int MAX_ITER = 10000;
double tolerance = 1e-15;
double pi = 4.0*atan(1.0);
int NI, NJ, ndims, niter;
double dt;
FILE *in;
char ip_file_type[50];
strcpy(ip_file_type,"ascii");
printf("Reading file \"solver.inp\"...\n");
in = fopen("solver.inp","r");
if (!in) {
fprintf(stderr,"Error: Input file \"solver.inp\" not found.\n");
return(0);
} 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);
fscanf(in,"%d",&NJ);
} 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 {
fprintf(stderr,"Error: Illegal format in solver.inp. Crash and burn!\n");
return(0);
}
}
fclose(in);
if (ndims != 2) {
fprintf(stderr,"ndims is not 2 in solver.inp. this code is to generate 2D initial conditions\n");
return(0);
}
printf("Grid: %d, %d\n",NI,NJ);
double tf = ((double)niter) * dt;
printf("Final Time: %lf\n",tf);
int i,j;
double dx = 1.0 / ((double)NI);
double dy = 1.0 / ((double)NJ);
{
double *x, *y, *u;
x = (double*) calloc (NI , sizeof(double));
y = (double*) calloc (NJ , sizeof(double));
u = (double*) calloc (NI*NJ, sizeof(double));
for (i = 0; i < NI; i++){
for (j = 0; j < NJ; j++){
x[i] = i*dx;
y[j] = j*dy;
int p = NJ*i + j;
u[p] = sin(2*pi*y[j]) * sin(2*pi*x[i]) / (2*ts*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,"%1.16e ",x[i]);
fprintf(out,"\n");
for (j = 0; j < NJ; j++) fprintf(out,"%1.16e ",y[j]);
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%1.16e ",u[p]);
}
}
fprintf(out,"\n");
fclose(out);
} else if ((!strcmp(ip_file_type,"binary")) || (!strcmp(ip_file_type,"bin"))) {
printf("Writing binary exact solution file initial.inp\n");
out = fopen("initial.inp","wb");
fwrite(x,sizeof(double),NI,out);
fwrite(y,sizeof(double),NJ,out);
double *U = (double*) calloc (NI*NJ,sizeof(double));
for (i=0; i < NI; i++) {
for (j = 0; j < NJ; j++) {
int p = NJ*i + j;
int q = NI*j + i;
U[q+0] = u[p];
}
}
fwrite(U,sizeof(double),NI*NJ,out);
free(U);
fclose(out);
}
free(x);
free(y);
free(u);
}
if (tf < ts) {
double *x, *y, *u;
x = (double*) calloc (NI , sizeof(double));
y = (double*) calloc (NJ , sizeof(double));
u = (double*) calloc (NI*NJ, sizeof(double));
for (i = 0; i < NI; i++){
for (j = 0; j < NJ; j++){
x[i] = i*dx;
y[j] = j*dy;
int p = NJ*i + j;
u[p] = sin(2*pi*y[j]) * sin(2*pi*x[i]) / (2*ts*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++){
for (j = 0; j < NJ; j++){
int p = NJ*i + j;
double new_u = sin(2.0*pi*(y[j]-u[p]*tf)) * sin(2.0*pi*(x[i]-u[p]*tf)) / (ts*2.0*pi);
double res = sqrt((new_u-u[p])*(new_u-u[p]));
if (res > maxres) maxres = res;
u[p] = new_u;
}
}
printf(" iter=%6d, max res=%1.6e\n", k, maxres);
if (maxres < tolerance) break;
}
FILE *out;
if (!strcmp(ip_file_type,"ascii")) {
out = fopen("exact.inp","w");
printf("Writing ASCII exact solution file exact.inp\n");
for (i = 0; i < NI; i++) fprintf(out,"%1.16e ",x[i]);
fprintf(out,"\n");
for (j = 0; j < NJ; j++) fprintf(out,"%1.16e ",y[j]);
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%1.16e ",u[p]);
}
}
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(y,sizeof(double),NJ,out);
double *U = (double*) calloc (NI*NJ,sizeof(double));
for (i=0; i < NI; i++) {
for (j = 0; j < NJ; j++) {
int p = NJ*i + j;
int q = NI*j + i;
U[q+0] = u[p];
}
}
fwrite(U,sizeof(double),NI*NJ,out);
free(U);
fclose(out);
}
free(x);
free(y);
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:
Note that iproc is set to
2 2
in solver.inp (i.e., 2 processors along x, and 2 processors along y). Thus, this example should be run with 4 MPI ranks (or change iproc).
After running the code, there should be 5 output files op_00000.dat, op_00001.dat, ... op_00004.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.
HyPar::op_file_format is set to tecplot2d in solver.inp, and thus, all the files are in a format that Tecplot (http://www.tecplot.com/) or other visualization software supporting the Tecplot format (e.g. VisIt - https://wci.llnl.gov/simulation/computer-codes/visit/) can read. In these files, the first two lines are the Tecplot headers, after which the data is written out as: the first two columns are grid indices, the next two columns are x and y coordinates, and the final column is the solution. HyPar::op_file_format can be set to text to get the solution files in plain text format (which can be read in and visualized in MATLAB for example).
The following figures show the solutions at each output time:
- t = 0 (initial)
- t = 0.4)
- t = 0.8
- t = 1.2
- t = 1.6 (final)
Since the exact solution is available at the final time the errors are calculated and reported on screen (see below) as well as errors.dat:
64 64 2 2 1.0000000000000001E-01 9.1090712454551764E-04 3.0493530313174199E-03 1.6532534290490630E-02 2.2016900000000000E-01 2.3717800000000000E-01
The numbers are: number of grid points in each dimension (HyPar::dim_global), number of processors in each dimension (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.
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:
64 64 2 2 1.0000000000000001E-01 1.5612974504934294E-17
The numbers are: number of grid points in each dimension (HyPar::dim_global), number of processors in each dimension (MPIVariables::iproc), time step size (HyPar::dt), and conservation error (HyPar::ConservationError).
Expected screen output:
HyPar - Parallel (MPI) version with 4 processes
Compiled with PETSc time integration.
Allocated simulation object(s).
Reading solver inputs from file "solver.inp".
No. of dimensions : 2
No. of variables : 1
Domain size : 64 64
Processes along each dimension : 2 2
No. of ghosts pts : 3
No. of iter. : 16
Restart iteration : 0
Time integration scheme : rk (44)
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 : 4
Initial solution file type : binary
Initial solution read mode : serial
Solution file write mode : serial
Solution file format : tecplot2d
Overwrite solution file : no
Physical model : burgers
Partitioning domain and allocating data arrays.
Reading array from binary file initial.inp (Serial mode).
Volume integral of the initial solution:
0: 1.7347234759768071E-18
Reading boundary conditions from boundary.inp.
Boundary periodic: Along dimension 0 and face +1
Boundary periodic: Along dimension 0 and face -1
Boundary periodic: Along dimension 1 and face +1
Boundary periodic: Along dimension 1 and face -1
4 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 16 iterations)
Writing solution file op_00000.dat.
Iteration: 1 Time: 1.000E-01 Max CFL: 5.093E-01 Max Diff. No.: -1.000E+00 Norm: 1.2181E-03 Conservation loss: 3.4694E-18
Iteration: 2 Time: 2.000E-01 Max CFL: 5.081E-01 Max Diff. No.: -1.000E+00 Norm: 1.2201E-03 Conservation loss: 8.6736E-18
Iteration: 3 Time: 3.000E-01 Max CFL: 5.093E-01 Max Diff. No.: -1.000E+00 Norm: 1.2242E-03 Conservation loss: 6.9389E-18
Iteration: 4 Time: 4.000E-01 Max CFL: 5.082E-01 Max Diff. No.: -1.000E+00 Norm: 1.2305E-03 Conservation loss: 1.2143E-17
Writing solution file op_00001.dat.
Iteration: 5 Time: 5.000E-01 Max CFL: 5.093E-01 Max Diff. No.: -1.000E+00 Norm: 1.2390E-03 Conservation loss: 8.6736E-18
Iteration: 6 Time: 6.000E-01 Max CFL: 5.083E-01 Max Diff. No.: -1.000E+00 Norm: 1.2500E-03 Conservation loss: 1.5612E-17
Iteration: 7 Time: 7.000E-01 Max CFL: 5.093E-01 Max Diff. No.: -1.000E+00 Norm: 1.2636E-03 Conservation loss: 1.7347E-18
Iteration: 8 Time: 8.000E-01 Max CFL: 5.083E-01 Max Diff. No.: -1.000E+00 Norm: 1.2803E-03 Conservation loss: 3.4694E-18
Writing solution file op_00002.dat.
Iteration: 9 Time: 9.000E-01 Max CFL: 5.093E-01 Max Diff. No.: -1.000E+00 Norm: 1.3003E-03 Conservation loss: 6.9388E-18
Iteration: 10 Time: 1.000E+00 Max CFL: 5.084E-01 Max Diff. No.: -1.000E+00 Norm: 1.3244E-03 Conservation loss: 6.9390E-18
Iteration: 11 Time: 1.100E+00 Max CFL: 5.093E-01 Max Diff. No.: -1.000E+00 Norm: 1.3532E-03 Conservation loss: 3.9899E-17
Iteration: 12 Time: 1.200E+00 Max CFL: 5.085E-01 Max Diff. No.: -1.000E+00 Norm: 1.3877E-03 Conservation loss: 1.3878E-17
Writing solution file op_00003.dat.
Iteration: 13 Time: 1.300E+00 Max CFL: 5.092E-01 Max Diff. No.: -1.000E+00 Norm: 1.4294E-03 Conservation loss: 1.3878E-17
Iteration: 14 Time: 1.400E+00 Max CFL: 5.086E-01 Max Diff. No.: -1.000E+00 Norm: 1.4802E-03 Conservation loss: 1.7347E-18
Iteration: 15 Time: 1.500E+00 Max CFL: 5.092E-01 Max Diff. No.: -1.000E+00 Norm: 1.5423E-03 Conservation loss: 4.6322E-22
Iteration: 16 Time: 1.600E+00 Max CFL: 5.086E-01 Max Diff. No.: -1.000E+00 Norm: 1.6190E-03 Conservation loss: 1.5613E-17
Writing solution file op_00004.dat.
Completed time integration (Final time: 1.600000).
Reading array from binary file exact.inp (Serial mode).
Computed errors for domain 0:
L1 Error : 9.1090712454551764E-04
L2 Error : 3.0493530313174199E-03
Linfinity Error : 1.6532534290490630E-02
Conservation Errors:
1.5612974504934294E-17
Solver runtime (in seconds): 2.2016900000000000E-01
Total runtime (in seconds): 2.3717800000000000E-01
Deallocating arrays.
Finished.
