HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
2D Linear Advection - Gaussian Pulse

Location: hypar/Examples/2D/LinearAdvection/GaussianPulse (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: 2D Linear Advection Equation (linearadr.h)

Domain: \(-6 \le x < 6, -3 \le y < 3\), "periodic" (_PERIODIC_) boundary conditions on all boundaries.

Initial solution: \(u\left(x,y,0\right) = u_0\left(x,y\right)= \exp\left[-\left(\frac{x^2}{2}+\frac{y^2}{2}\right)\right]\)
Exact solution: \(u\left(x,y,t\right) = u_0\left(x-a_xt,y-a_yt\right)\).

Numerical Method:

Input files required:

solver.inp

begin
ndims 2
nvars 1
size 120 60
iproc 4 2
ghost 3
n_iter 300
time_scheme rk
time_scheme_type ssprk3
hyp_space_scheme crweno5
conservation_check yes
dt 0.04
screen_op_iter 5
file_op_iter 15
op_file_format tecplot2d
op_overwrite no
model linear-advection-diffusion-reaction
end

boundary.inp

4
periodic 0 1 0 0 -3.0 3.0
periodic 0 -1 0 0 -3.0 3.0
periodic 1 1 -6.0 6.0 0 0
periodic 1 -1 -6.0 6.0 0 0

physics.inp (specifies \(a_x\) and \(a_y\))

begin
advection 1.0 0
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, compile and run the following code in the run directory. Note: if the final time is an integer multiple of the time period, the file initial.inp can also be used as the exact solution exact.inp (i.e. create a sym link called exact.inp pointing to initial.inp, or just copy initial.inp to exact.inp).

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
int main(){
int NI=100,NJ=50,ndims=2;
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);
fscanf(in,"%d",&NJ);
} 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 != 2) {
printf("ndims is not 2 in solver.inp. this code is to generate 2D initial conditions\n");
return(0);
}
printf("Grid:\t\t\t%d X %d\n",NI,NJ);
int i,j;
double dx = 12.0 / ((double)NI);
double dy = 6.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] = -6 + i*dx;
y[j] = -3 + j*dy;
int p = NJ*i + j;
u[p] = exp(-((x[i]*x[i]/2+y[j]*y[j]/2)));
}
}
FILE *out;
if (!strcmp(ip_file_type,"ascii")) {
out = fopen("initial.inp","w");
for (i = 0; i < NI; i++) fprintf(out,"%lf ",x[i]);
fprintf(out,"\n");
for (j = 0; j < NJ; j++) fprintf(out,"%lf ",y[j]);
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%lf ",u[p]);
}
}
fprintf(out,"\n");
fclose(out);
} else if ((!strcmp(ip_file_type,"binary")) || (!strcmp(ip_file_type,"bin"))) {
printf("Error: Writing binary initial solution file not implemented. ");
printf("Please choose ip_file_type in solver.inp as \"ascii\".\n");
}
free(x);
free(y);
free(u);
return(0);
}

Output:

Note that iproc is set to 

  4 2

in solver.inp (i.e., 4 processors along x, and 2 processors along y). Thus, this example should be run with 8 MPI ranks (or change iproc).

After running the code, there should be 21 output files op_00000.dat, op_00001.dat, ... op_00020.dat; the first one is the solution at \(t=0\) and the final one is the solution at \(t=12\). 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 animation was generated from the solution files:

Solution_2DLinearAdvGauss.gif

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

120 60 4 2 4.0000000000000001E-02 9.2322365120429538E-05 8.4380658026143653E-05 9.9005102091931363E-05 2.3760400000000002E+00 2.3815119999999999E+00

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:

120 60 4 2 4.0000000000000001E-02 4.2523145153349166E-16

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 8 processes
Compiled with PETSc time integration.
Reading solver inputs from file "solver.inp".
No. of dimensions : 2
No. of variables : 1
Domain size : 120 60
Processes along each dimension : 4 2
No. of ghosts pts : 3
No. of iter. : 300
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 : 4.000000E-02
Check for conservation : yes
Screen output iterations : 5
File output iterations : 15
Initial solution file type : ascii
Initial solution read mode : serial
Solution file write mode : serial
Solution file format : tecplot2d
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: 6.2660822699999983E+00
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 = "linear-advection-diffusion-reaction"
Reading physical model inputs from file "physics.inp".
Setting up time integration.
Solving in time (from 0 to 300 iterations)
Writing solution file op_00000.dat.
Iteration: 5 Time: 2.000E-01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9075E-03 Conservation loss: 2.8349E-16
Iteration: 10 Time: 4.000E-01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9075E-03 Conservation loss: 8.5046E-16
Iteration: 15 Time: 6.000E-01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9075E-03 Conservation loss: 2.8349E-16
Writing solution file op_00001.dat.
Iteration: 20 Time: 8.000E-01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9075E-03 Conservation loss: 1.4174E-16
Iteration: 25 Time: 1.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9075E-03 Conservation loss: 5.6698E-16
Iteration: 30 Time: 1.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9075E-03 Conservation loss: 5.6698E-16
Writing solution file op_00002.dat.
Iteration: 35 Time: 1.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 1.4174E-16
Iteration: 40 Time: 1.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 2.8349E-16
Iteration: 45 Time: 1.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 2.8349E-16
Writing solution file op_00003.dat.
Iteration: 50 Time: 2.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 0.0000E+00
Iteration: 55 Time: 2.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 2.8349E-16
Iteration: 60 Time: 2.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 8.5046E-16
Writing solution file op_00004.dat.
Iteration: 65 Time: 2.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 1.4174E-16
Iteration: 70 Time: 2.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9074E-03 Conservation loss: 2.8349E-16
Iteration: 75 Time: 3.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 4.2523E-16
Writing solution file op_00005.dat.
Iteration: 80 Time: 3.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 4.2523E-16
Iteration: 85 Time: 3.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 1.4174E-16
Iteration: 90 Time: 3.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 5.6698E-16
Writing solution file op_00006.dat.
Iteration: 95 Time: 3.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 7.0872E-16
Iteration: 100 Time: 4.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 1.4174E-16
Iteration: 105 Time: 4.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 0.0000E+00
Writing solution file op_00007.dat.
Iteration: 110 Time: 4.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9073E-03 Conservation loss: 5.6698E-16
Iteration: 115 Time: 4.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 1.4174E-16
Iteration: 120 Time: 4.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 5.6698E-16
Writing solution file op_00008.dat.
Iteration: 125 Time: 5.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 5.6698E-16
Iteration: 130 Time: 5.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 2.8349E-16
Iteration: 135 Time: 5.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 1.4174E-16
Writing solution file op_00009.dat.
Iteration: 140 Time: 5.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 1.2757E-15
Iteration: 145 Time: 5.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 4.2523E-16
Iteration: 150 Time: 6.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9072E-03 Conservation loss: 0.0000E+00
Writing solution file op_00010.dat.
Iteration: 155 Time: 6.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 4.2523E-16
Iteration: 160 Time: 6.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 1.1340E-15
Iteration: 165 Time: 6.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 1.1340E-15
Writing solution file op_00011.dat.
Iteration: 170 Time: 6.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 2.8349E-16
Iteration: 175 Time: 7.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 1.4174E-16
Iteration: 180 Time: 7.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 5.6698E-16
Writing solution file op_00012.dat.
Iteration: 185 Time: 7.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 5.6698E-16
Iteration: 190 Time: 7.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 1.4174E-15
Iteration: 195 Time: 7.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9071E-03 Conservation loss: 7.0872E-16
Writing solution file op_00013.dat.
Iteration: 200 Time: 8.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 5.6698E-16
Iteration: 205 Time: 8.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 4.2523E-16
Iteration: 210 Time: 8.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 2.8349E-16
Writing solution file op_00014.dat.
Iteration: 215 Time: 8.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 4.2523E-16
Iteration: 220 Time: 8.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 0.0000E+00
Iteration: 225 Time: 9.000E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 1.4174E-16
Writing solution file op_00015.dat.
Iteration: 230 Time: 9.200E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 1.4174E-16
Iteration: 235 Time: 9.400E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9070E-03 Conservation loss: 7.0872E-16
Iteration: 240 Time: 9.600E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 1.4174E-16
Writing solution file op_00016.dat.
Iteration: 245 Time: 9.800E+00 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 1.4174E-16
Iteration: 250 Time: 1.000E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 1.4174E-16
Iteration: 255 Time: 1.020E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 4.2523E-16
Writing solution file op_00017.dat.
Iteration: 260 Time: 1.040E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 2.8349E-16
Iteration: 265 Time: 1.060E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 4.2523E-16
Iteration: 270 Time: 1.080E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 1.4174E-16
Writing solution file op_00018.dat.
Iteration: 275 Time: 1.100E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9069E-03 Conservation loss: 8.5046E-16
Iteration: 280 Time: 1.120E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9068E-03 Conservation loss: 1.4174E-16
Iteration: 285 Time: 1.140E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9068E-03 Conservation loss: 7.0872E-16
Writing solution file op_00019.dat.
Iteration: 290 Time: 1.160E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9068E-03 Conservation loss: 2.8349E-16
Iteration: 295 Time: 1.180E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9068E-03 Conservation loss: 0.0000E+00
Iteration: 300 Time: 1.200E+01 Max CFL: 4.000E-01 Max Diff. No.: 0.000E+00 Norm: 5.9068E-03 Conservation loss: 4.2523E-16
Writing solution file op_00020.dat.
Completed time integration (Final time: 12.000000).
Reading array from ASCII file exact.inp (Serial mode).
Computed errors:
L1 Error : 9.2322365120429538E-05
L2 Error : 8.4380658026143653E-05
Linfinity Error : 9.9005102091931363E-05
Conservation Errors:
4.2523145153349166E-16
Solver runtime (in seconds): 2.3760400000000002E+00
Total runtime (in seconds): 2.3815119999999999E+00
Deallocating arrays.
Finished.