See 2D Euler Equations - Isentropic Vortex Convection to familiarize yourself with this case. This example uses a DMD object that has already been trained (see 2D Euler Equations - Isentropic Vortex Convection).
Location: hypar/Examples/2D/NavierStokes2D/InviscidVortexConvection_libROM_DMD_Predict (This directory contains all the input files needed to run this case.)
Reduced Order Modeling: This example predicts the solution from a trained DMD object. The code does not solve the PDE by discretizing in space and integrating in time.
Governing equations: 2D Euler Equations (navierstokes2d.h - By default, NavierStokes2D::Re is set to -1 which makes the code skip the parabolic terms, i.e., the 2D Euler equations are solved.)
Reference: C.-W. Shu, "Essentially Non-oscillatory and Weighted Essentially
Non-oscillatory Schemes for Hyperbolic Conservation Laws", ICASE Report 97-65, 1997
Domain: \(0 \le x,y \le 10\), "periodic" (_PERIODIC_) boundary conditions.
Initial solution: The freestream flow is given by
\begin{equation} \rho_\infty = 1,\ u_\infty = 0.1,\ v_\infty = 0,\ p_\infty = 1 \end{equation}
and a vortex is introduced, specified as
\begin{align} \rho &= \left[ 1 - \frac{\left(\gamma-1\right)b^2}{8\gamma\pi^2} e^{1-r^2} \right]^{\frac{1}{\gamma-1}},\ p = \rho^\gamma, \\ u &= u_\infty - \frac{b}{2\pi} e^{\frac{1}{2}\left(1-r^2\right)} \left(y-y_c\right),\ v = v_\infty + \frac{b}{2\pi} e^{\frac{1}{2}\left(1-r^2\right)} \left(x-x_c\right), \end{align}
where \(b=0.5\) is the vortex strength and \(r = \left[(x-x_c)^2 + (y-y_c)^2 \right]^{1/2}\) is the distance from the vortex center \(\left(x_c,y_c\right) = \left(5,5\right)\).
Reduced Order Modeling:
Note:
In this mode, HyPar will run just like an usual PDE simulation, except that it will swap out the numerical spatial discretization and time integration with the ROM-based prediction. The input files and output files will be the same as a regular simulation with the following comments:
- Numerical method inputs are ignored (eg. those that specify spatial discretization and time integration methods).
- In solver.inp, the values for dt, n_iter, and file_op_iter is used only to compute the simulation times at which to compute and write the solution. The time step size, dt, need not respect any CFL criterion.
- HyPar::ConservationCheck is set to "no" since it is not possible to compute conservation loss for a general domain (because boundary fluxes are not being computed).
Input files required:
librom.inp
begin
mode predict
dmd_dirname DMD
end
DMD Object(s) :
The trained DMD object(s) must be located in the directory specified in librom.inp as dmd_dirname (DMDROMObject::m_dirname). For this example, they were generated using 2D Euler Equations - Isentropic Vortex Convection.
solver.inp
begin
ndims 2
nvars 4
size 64 64
iproc 4 4
ghost 3
n_iter 40
restart_iter 0
dt 1.0
file_op_iter 1
ip_file_type binary
op_file_format binary
op_overwrite no
model navierstokes2d
end
boundary.inp
4
periodic 0 1 0 0 0 10.0
periodic 0 -1 0 0 0 10.0
periodic 1 1 0 10.0 0 0
periodic 1 -1 0 10.0 0 0
physics.inp
begin
gamma 1.4
upwinding roe
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 power(double x,double a)
{
return(exp(a*log(x)));
}
const double pi = 4.0*atan(1.0);
const double GAMMA = 1.4;
int NI,NJ,ndims,n_iter;
double tf, dt;
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");
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, "n_iter")) fscanf(in,"%d",&n_iter);
else if (!strcmp(word, "dt")) fscanf(in,"%lf",&dt);
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 exact solution\n");
return(0);
}
printf("Grid:\t\t\t%d X %d\n",NI,NJ);
int i,j;
double dx = 10.0 / ((double)NI);
double dy = 10.0 / ((double)NJ);
tf = (double)n_iter * dt; double tff = tf;
while (tf > 20) tf -= 20;
double *x, *y, *u0, *u1, *u2, *u3;
x = (double*) calloc (NI , sizeof(double));
y = (double*) calloc (NJ , sizeof(double));
u0 = (double*) calloc (NI*NJ, sizeof(double));
u1 = (double*) calloc (NI*NJ, sizeof(double));
u2 = (double*) calloc (NI*NJ, sizeof(double));
u3 = (double*) calloc (NI*NJ, sizeof(double));
double u_inf = 0.5;
double v_inf = 0.0;
double b = u_inf;
double x0, y0;
x0 = 5.0, y0 = 5.0;
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;
double rx, ry;
rx = (x[i] - x0);
ry = (y[j] - y0);
if (rx < -5) { rx += 10; }
else if (rx > 5) { rx -= 10; }
double rsq = rx*rx + ry*ry;
double rho, u, v, P;
double du, dv;
rho = power(1.0 - ((GAMMA-1.0)*b*b)/(8.0*GAMMA*pi*pi) * exp(1.0-rsq), 1.0/(GAMMA-1.0));
P = power(rho,GAMMA);
du = - b/(2.0*pi) * exp(0.5*(1.0-rsq)) * ry;
dv = b/(2.0*pi) * exp(0.5*(1.0-rsq)) * rx;
u = u_inf + du;
v = v_inf + dv;
u0[p] = rho;
u1[p] = rho*u;
u2[p] = rho*v;
u3[p] = P/(GAMMA-1.0) + 0.5*rho*(u*u+v*v);
}
}
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 (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 ",u0[p]);
}
}
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%lf ",u1[p]);
}
}
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%lf ",u2[p]);
}
}
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%lf ",u3[p]);
}
}
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(y,sizeof(double),NJ,out);
double *U = (double*) calloc (4*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[4*q+0] = u0[p];
U[4*q+1] = u1[p];
U[4*q+2] = u2[p];
U[4*q+3] = u3[p];
}
}
fwrite(U,sizeof(double),4*NI*NJ,out);
free(U);
fclose(out);
}
x0 = 5.0+tf*u_inf, y0 = 5.0;
if (x0 > 10) x0 -= 10;
printf("Final time: %lf, Vortex center: %lf, %lf\n",tff,x0,y0);
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;
double rx, ry;
rx = (x[i] - x0);
ry = (y[j] - y0);
if (rx < -5) { rx += 10; }
else if (rx > 5) { rx -= 10; }
double rsq = rx*rx + ry*ry;
double rho, u, v, P;
double du, dv;
rho = power(1.0 - ((GAMMA-1.0)*b*b)/(8.0*GAMMA*pi*pi) * exp(1.0-rsq), 1.0/(GAMMA-1.0));
P = power(rho,GAMMA);
du = - b/(2.0*pi) * exp(0.5*(1.0-rsq)) * ry;
dv = b/(2.0*pi) * exp(0.5*(1.0-rsq)) * rx;
u = u_inf + du;
v = v_inf + dv;
u0[p] = rho;
u1[p] = rho*u;
u2[p] = rho*v;
u3[p] = P/(GAMMA-1.0) + 0.5*rho*(u*u+v*v);
}
}
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 (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 ",u0[p]);
}
}
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%lf ",u1[p]);
}
}
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%lf ",u2[p]);
}
}
fprintf(out,"\n");
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = NJ*i + j;
fprintf(out,"%lf ",u3[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 (4*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[4*q+0] = u0[p];
U[4*q+1] = u1[p];
U[4*q+2] = u2[p];
U[4*q+3] = u3[p];
}
}
fwrite(U,sizeof(double),4*NI*NJ,out);
free(U);
fclose(out);
}
free(x);
free(y);
free(u0);
free(u1);
free(u2);
free(u3);
return(0);
}
Output:
Note that iproc is set to
4 4
in solver.inp (i.e., 4 processors along x, and 4 processors along y). Thus, this example should be run with 16 MPI ranks (or change iproc).
After running the code, there should be the following output files:
- 41 output files op_00000.bin, op_00001.bin, ... op_00040.bin; these are the solutions as predicted by the ROM .
The first of each of these file sets is the solution at \(t=0\) and the final one is the solution at \(t=40\). Since HyPar::op_overwrite is set to no in solver.inp, a separate file is written for solutions at each output time. All the files are binary (HyPar::op_file_format is set to binary in solver.inp).
The provided Python script (Examples/Python/plotSolution_2DBinary.py) can be used to generate plots from these binary files. Alternatively, HyPar::op_file_format can be set to text or tecplot2d, and Tecplot/VisIt or something similar can be used to plot the resulting files.
The animation shows the evolution of the solution (density).
Since the exact solution is available at the final time (exact.inp), the numerical errors are calculated and reported on screen (see below) as well as errors.dat:
64 64 4 4 1.0000000000000000E+00 1.7057960004704812E-06 4.1746751056680860E-06 2.4053620904213429E-05 3.8236610000000000E+00 4.0025190000000004E+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.
Expected screen output:
HyPar - Parallel (MPI) version with 16 processes
Compiled with PETSc time integration.
Allocated simulation object(s).
Reading solver inputs from file "solver.inp".
No. of dimensions : 2
No. of variables : 4
Domain size : 64 64
Processes along each dimension : 4 4
Exact solution domain size : 64 64
No. of ghosts pts : 3
No. of iter. : 40
Restart iteration : 0
Time integration scheme : euler
Spatial discretization scheme (hyperbolic) : 1
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+00
Check for conservation : no
Screen output iterations : 1
File output iterations : 1
Initial solution file type : binary
Initial solution read mode : serial
Solution file write mode : serial
Solution file format : binary
Overwrite solution file : no
Physical model : navierstokes2d
Partitioning domain and allocating data arrays.
Reading array from binary file initial.inp (Serial mode).
Volume integral of the initial solution:
0: 9.9980704056028685E+01
1: 4.9990352985508352E+01
2: -9.5749398777884853E-07
3: 2.6245709189715654E+02
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.
Initializing physics. Model = "navierstokes2d"
Reading physical model inputs from file "physics.inp".
Setting up libROM interface.
libROM inputs and parameters:
reduced model dimensionality: 0
sampling frequency: 0
mode: predict
type: DMD
save to file: true
libROM DMD inputs:
number of samples per window: 2147483647
directory name for DMD onjects: DMD
libROMInterface::loadROM() - loading ROM objects.
Loading DMD object (DMD/dmdobj_0000), time window=[0.00e+00,-1.00e+00].
libROM: Predicted solution at time 0.0000e+00 using ROM, wallclock time: 0.013085.
Writing solution file op_00000.bin.
libROM: Predicted solution at time 1.0000e+00 using ROM, wallclock time: 0.035576.
Writing solution file op_00001.bin.
libROM: Predicted solution at time 2.0000e+00 using ROM, wallclock time: 0.033418.
Writing solution file op_00002.bin.
libROM: Predicted solution at time 3.0000e+00 using ROM, wallclock time: 0.135173.
Writing solution file op_00003.bin.
libROM: Predicted solution at time 4.0000e+00 using ROM, wallclock time: 0.089607.
Writing solution file op_00004.bin.
libROM: Predicted solution at time 5.0000e+00 using ROM, wallclock time: 0.107563.
Writing solution file op_00005.bin.
libROM: Predicted solution at time 6.0000e+00 using ROM, wallclock time: 0.060770.
Writing solution file op_00006.bin.
libROM: Predicted solution at time 7.0000e+00 using ROM, wallclock time: 0.125400.
Writing solution file op_00007.bin.
libROM: Predicted solution at time 8.0000e+00 using ROM, wallclock time: 0.116264.
Writing solution file op_00008.bin.
libROM: Predicted solution at time 9.0000e+00 using ROM, wallclock time: 0.098884.
Writing solution file op_00009.bin.
libROM: Predicted solution at time 1.0000e+01 using ROM, wallclock time: 0.057265.
Writing solution file op_00010.bin.
libROM: Predicted solution at time 1.1000e+01 using ROM, wallclock time: 0.050350.
Writing solution file op_00011.bin.
libROM: Predicted solution at time 1.2000e+01 using ROM, wallclock time: 0.032244.
Writing solution file op_00012.bin.
libROM: Predicted solution at time 1.3000e+01 using ROM, wallclock time: 0.034619.
Writing solution file op_00013.bin.
libROM: Predicted solution at time 1.4000e+01 using ROM, wallclock time: 0.086340.
Writing solution file op_00014.bin.
libROM: Predicted solution at time 1.5000e+01 using ROM, wallclock time: 0.031885.
Writing solution file op_00015.bin.
libROM: Predicted solution at time 1.6000e+01 using ROM, wallclock time: 0.041468.
Writing solution file op_00016.bin.
libROM: Predicted solution at time 1.7000e+01 using ROM, wallclock time: 0.036360.
Writing solution file op_00017.bin.
libROM: Predicted solution at time 1.8000e+01 using ROM, wallclock time: 0.077352.
Writing solution file op_00018.bin.
libROM: Predicted solution at time 1.9000e+01 using ROM, wallclock time: 0.044662.
Writing solution file op_00019.bin.
libROM: Predicted solution at time 2.0000e+01 using ROM, wallclock time: 0.042875.
Writing solution file op_00020.bin.
libROM: Predicted solution at time 2.1000e+01 using ROM, wallclock time: 0.047680.
Writing solution file op_00021.bin.
libROM: Predicted solution at time 2.2000e+01 using ROM, wallclock time: 0.067208.
Writing solution file op_00022.bin.
libROM: Predicted solution at time 2.3000e+01 using ROM, wallclock time: 0.030006.
Writing solution file op_00023.bin.
libROM: Predicted solution at time 2.4000e+01 using ROM, wallclock time: 0.100523.
Writing solution file op_00024.bin.
libROM: Predicted solution at time 2.5000e+01 using ROM, wallclock time: 0.067925.
Writing solution file op_00025.bin.
libROM: Predicted solution at time 2.6000e+01 using ROM, wallclock time: 0.033365.
Writing solution file op_00026.bin.
libROM: Predicted solution at time 2.7000e+01 using ROM, wallclock time: 0.054427.
Writing solution file op_00027.bin.
libROM: Predicted solution at time 2.8000e+01 using ROM, wallclock time: 0.129765.
Writing solution file op_00028.bin.
libROM: Predicted solution at time 2.9000e+01 using ROM, wallclock time: 0.063883.
Writing solution file op_00029.bin.
libROM: Predicted solution at time 3.0000e+01 using ROM, wallclock time: 0.054281.
Writing solution file op_00030.bin.
libROM: Predicted solution at time 3.1000e+01 using ROM, wallclock time: 0.071597.
Writing solution file op_00031.bin.
libROM: Predicted solution at time 3.2000e+01 using ROM, wallclock time: 0.070207.
Writing solution file op_00032.bin.
libROM: Predicted solution at time 3.3000e+01 using ROM, wallclock time: 0.073919.
Writing solution file op_00033.bin.
libROM: Predicted solution at time 3.4000e+01 using ROM, wallclock time: 0.077781.
Writing solution file op_00034.bin.
libROM: Predicted solution at time 3.5000e+01 using ROM, wallclock time: 0.067210.
Writing solution file op_00035.bin.
libROM: Predicted solution at time 3.6000e+01 using ROM, wallclock time: 0.080262.
Writing solution file op_00036.bin.
libROM: Predicted solution at time 3.7000e+01 using ROM, wallclock time: 0.046820.
Writing solution file op_00037.bin.
libROM: Predicted solution at time 3.8000e+01 using ROM, wallclock time: 0.059902.
Writing solution file op_00038.bin.
libROM: Predicted solution at time 3.9000e+01 using ROM, wallclock time: 0.044101.
Writing solution file op_00039.bin.
libROM: Predicted solution at time 4.0000e+01 using ROM, wallclock time: 0.095935.
Writing solution file op_00040.bin.
Reading array from binary file exact.inp (Serial mode).
libROM: total prediction/query wallclock time: 2.687957 (seconds).
Computed errors for domain 0:
L1 Error : 1.7057960004704812E-06
L2 Error : 4.1746751056680860E-06
Linfinity Error : 2.4053620904213429E-05
Solver runtime (in seconds): 3.8236610000000000E+00
Total runtime (in seconds): 4.0025190000000004E+00
Deallocating arrays.
Finished.
