HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
3D Navier-Stokes - Inviscid Shock-Cylinder Interaction (Time-Windowed DMD)

See Inviscid Shock-Cylinder Interaction to familiarize yourself with this case.

Location: hypar/Examples/3D/NavierStokes3D/2D_Shock_Cylinder_Interaction_libROM_DMD_Train

Governing equations: 3D Navier-Stokes Equations (navierstokes3d.h - by default NavierStokes3D::Re is set to -1 which makes the code skip the parabolic terms, i.e., the 3D Euler equations are solved.)

Reduced Order Modeling: This example trains a DMD object and then predicts the solution using the DMD at the same times that the actual HyPar solution is written at.

Domain: \(-2.5 \le x \le 7.5\), \(-5 \le y \le 5\) Note: This is a 2D flow simulated using a 3-dimensional setup by taking the length of the domain along z to be very small and with only 3 grid points (the domain size along z must be smaller than the cylinder length).

Geometry: A cylinder of radius 1.0 centered at (0,0) (hypar/Examples/STLGeometries/cylinder.stl)

Boundary conditions:

Reference:

  • O. Boiron, G. Chiavassa, R. Donat, A high-resolution penalization method for large Mach number flows in the presence of obstacles, Computers & Fluids, 38 (2009), pp. 703-714, Section 5.1

Initial solution: see reference

Other parameters:

Numerical Method:

Reduced Order Modeling:

Input files required:

These files are all located in: hypar/Examples/3D/NavierStokes3D/2D_Shock_Cylinder_Interaction_libROM_DMD_Train/

librom.inp

begin
rdim 32
sampling_frequency 1
mode train
dmd_num_win_samples 100
end

solver.inp

begin
ndims 3
nvars 5
size 256 256 3
iproc 4 4 1
ghost 3
n_iter 400
time_scheme rk
time_scheme_type 44
hyp_space_scheme weno5
hyp_interp_type characteristic
dt 0.005
screen_op_iter 10
file_op_iter 40
ip_file_type binary
op_file_format binary
op_overwrite no
model navierstokes3d
immersed_body cylinder.stl
end

boundary.inp

6
subsonic-inflow 0 1 0 0 -9999.0 9999.0 -9999.0 9999.0
3.8571428571428572 2.6293687924887181 0.0 0.0
supersonic-inflow 0 -1 0 0 -9999.0 9999.0 -9999.0 9999.0
1.0 0.0 0.0 0.0 1.0
slip-wall 1 1 -9999.0 9999.0 0 0 -9999.0 9999.0
slip-wall 1 -1 -9999.0 9999.0 0 0 -9999.0 9999.0
periodic 2 1 -9999.0 9999.0 -9999.0 9999.0 0 0
periodic 2 -1 -9999.0 9999.0 -9999.0 9999.0 0 0

physics.inp

begin
gamma 1.4
upwinding roe
end

cylinder.stl : the filename "cylinder.stl" must match the input for immersed_body in solver.inp.
Located at hypar/Examples/STLGeometries/cylinder.stl

To generate initial.inp (initial 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()
{
const double gamma = 1.4;
int NI,NJ,NK,ndims;
char ip_file_type[50];
FILE *in;
strcpy(ip_file_type,"ascii");
printf("Reading file \"solver.inp\"...\n");
in = fopen("solver.inp","r");
if (!in) {
printf("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);
fscanf(in,"%d",&NK);
} 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");
return(0);
}
}
fclose(in);
if (ndims != 3) {
printf("Error: ndims is not 3 in solver.inp.\n");
return(0);
}
printf("Grid: %3d x %3d x %3d.\n",NI,NJ,NK);
double Lx = 10.0; /* length of domain along x */
double Ly = 10.0; /* length of domain along y */
double xmin = -2.5; /* left x-boundary */
double ymin = -5.0; /* left y-boundary */
double xs = -2.0; /* location of initial shock */
double Ms = 3.0; /* shock Mach number */
int i,j,k;
double *x, *y, *z;
double dx, dy, dz;
x = (double*) calloc (NI, sizeof(double));
y = (double*) calloc (NJ, sizeof(double));
z = (double*) calloc (NK, sizeof(double));
dx = Lx / ((double)(NI-1));
for (i=0; i<NI; i++) x[i] = xmin + i*dx;
dy = Ly / ((double)(NJ-1));
for (j=0; j<NJ; j++) y[j] = ymin + j*dy;
dz = 0.5 * (dx +dy);
for (k = 0; k < NK; k++) z[k] = (k-NK/2) * dz;
/* pre-shock conditions */
double rho_inf = 1.0;
double p_inf = 1.0;
double u_inf = 0.0; //- Ms * sqrt(gamma*p_inf/rho_inf);
double v_inf = 0.0;
double w_inf = 0.0;
/* post-shock conditions */
double rho_ps = rho_inf * ((gamma+1.0)*Ms*Ms) / ((gamma-1.0)*Ms*Ms+2.0);
double p_ps = p_inf * (2.0*gamma*Ms*Ms-(gamma-1.0)) / (gamma+1.0);
double M2 = sqrt (((gamma-1.0)*Ms*Ms+2.0) / (2.0*gamma*Ms*Ms-(gamma-1.0)));
double u_ps = - (M2 * sqrt(gamma*p_ps/rho_ps) - Ms * sqrt(gamma*p_inf/rho_inf));
double v_ps = 0.0;
double w_ps = 0.0;
printf("Pre-shock conditions:-\n");
printf("Density = %1.16E\n", rho_inf);
printf("Pressure = %1.16E\n", p_inf );
printf("u-velocity = %1.16E\n", u_inf );
printf("v-velocity = %1.16E\n", v_inf );
printf("w-velocity = %1.16E\n", w_inf );
printf("Mach number= %1.16E\n", Ms );
printf("Post-shock conditions:-\n");
printf("Density = %1.16E\n", rho_ps);
printf("Pressure = %1.16E\n", p_ps );
printf("u-velocity = %1.16E\n", u_ps );
printf("v-velocity = %1.16E\n", v_ps );
printf("w-velocity = %1.16E\n", w_ps );
printf("Mach number= %1.16E\n", M2 );
printf("Ratios:-\n");
printf("p2/p1 = %1.16E\n",p_ps/p_inf);
printf("rho2/rho1 = %1.16E\n",rho_ps/rho_inf);
double *U = (double*) calloc (5*NI*NJ*NK, sizeof(double));
for (i = 0; i < NI; i++){
for (j = 0; j < NJ; j++){
for (k = 0; k < NK; k++){
int p = i + NI*j + NI*NJ*k;
double rho, u, v, w, P;
if (x[i] > xs) {
rho = rho_inf;
u = u_inf;
v = v_inf;
w = w_inf;
P = p_inf;
} else {
rho = rho_ps;
u = u_ps;
v = v_ps;
w = w_ps;
P = p_ps;
}
U[5*p+0] = rho;
U[5*p+1] = rho*u;
U[5*p+2] = rho*v;
U[5*p+3] = rho*w;
U[5*p+4] = P/(gamma-1.0) + 0.5*rho*(u*u+v*v+w*w);
}
}
}
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 (k = 0; k < NK; k++) fprintf(out,"%1.16E ",z[k]);
fprintf(out,"\n");
for (k = 0; k < NK; k++) {
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = i + NK*j + NI*NJ*k;
fprintf(out,"%1.16E ",U[5*p+0]);
}
}
}
fprintf(out,"\n");
for (k = 0; k < NK; k++) {
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = i + NK*j + NI*NJ*k;
fprintf(out,"%1.16E ",U[5*p+1]);
}
}
}
fprintf(out,"\n");
for (k = 0; k < NK; k++) {
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = i + NK*j + NI*NJ*k;
fprintf(out,"%1.16E ",U[5*p+2]);
}
}
}
fprintf(out,"\n");
for (k = 0; k < NK; k++) {
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = i + NK*j + NI*NJ*k;
fprintf(out,"%1.16E ",U[5*p+3]);
}
}
}
fprintf(out,"\n");
for (k = 0; k < NK; k++) {
for (j = 0; j < NJ; j++) {
for (i = 0; i < NI; i++) {
int p = i + NK*j + NI*NJ*k;
fprintf(out,"%1.16E ",U[5*p+4]);
}
}
}
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);
fwrite(z,sizeof(double),NK,out);
fwrite(U,sizeof(double),5*NI*NJ*NK,out);
fclose(out);
}
free(x);
free(y);
free(z);
free(U);
return(0);
}

Output:

Note that iproc is set to

  4 4 1

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

Please see the original example in the "Immersed Boundaries Examples" section for a full description of the output files and how to view them. The following only contains libROM-specific comments.

After running the code, there should be the following output files:

  • 11 output file op_00000.bin, ..., op_00010.bin; these are the HyPar solutions.
  • 11 output file op_rom_00000.bin, ..., op_rom_00010.bin; these are the predicted solutions from the DMD object(s).

All the files are binary (HyPar::op_file_format is set to binary in solver.inp).

The provided Python script (plotSolution.py) can be used to generate plots from the binary files that compare the HyPar and DMD solutions. Alternatively, HyPar::op_file_format can be set to tecplot3d, and Tecplot/VisIt or something similar can be used to plot the resulting text files.

The following plot shows the density contours for the final solution (t=2). FOM (full-order model) refers to the HyPar solution, ROM (reduced-order model) refers to the DMD solution, and Diff is the difference between the two.

Solution_3DNavStokShockCyl_Density_libROM_DMD.png

Wall clock times:

  • PDE solution: 815 seconds
  • DMD training time: 155 seconds
  • DMD prediction/query time: 33 seconds

The L1, L2, and Linf norms of the diff between the HyPar and ROM solution at the final time are calculated and reported on screen (see below) as well as pde_rom_diff.dat:

256 256 3 4 4 1 5.0000000000000001E-03 6.4931422313752370E-03 1.9356397835829150E-02 9.1818439088445988E-02 1.0099263970000000E+03 1.0104446960000000E+03

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 norms of the diff (HyPar::rom_diff_norms), 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 : 3
No. of variables : 5
Domain size : 256 256 3
Processes along each dimension : 4 4 1
Exact solution domain size : 256 256 3
No. of ghosts pts : 3
No. of iter. : 400
Restart iteration : 0
Time integration scheme : rk (44)
Spatial discretization scheme (hyperbolic) : weno5
Split hyperbolic flux term? : no
Interpolation type for hyperbolic term : characteristic
Spatial discretization type (parabolic ) : nonconservative-1stage
Spatial discretization scheme (parabolic ) : 2
Time Step : 5.000000E-03
Check for conservation : no
Screen output iterations : 10
File output iterations : 40
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 : navierstokes3d
Immersed Body : cylinder.stl
Partitioning domain and allocating data arrays.
Reading array from binary file initial.inp (Serial mode).
Volume integral of the initial solution:
0: 1.3577505742780760E+01
1: 6.1066251110367062E+00
2: 0.0000000000000000E+00
3: 0.0000000000000000E+00
4: 5.1720680582882707E+01
Reading boundary conditions from boundary.inp.
Boundary subsonic-inflow: Along dimension 0 and face +1
Boundary supersonic-inflow: Along dimension 0 and face -1
Boundary slip-wall: Along dimension 1 and face +1
Boundary slip-wall: Along dimension 1 and face -1
Boundary periodic: Along dimension 2 and face +1
Boundary periodic: Along dimension 2 and face -1
6 boundary condition(s) read.
Immersed body read from cylinder.stl:
Number of facets: 628
Bounding box: [-1.0,1.0] X [-1.0,1.0] X [-5.0,5.0]
Number of grid points inside immersed body: 6408 ( 3.3%).
Number of immersed boundary points : 1242 ( 0.6%).
Immersed body simulation mode : 2d (xy).
Initializing solvers.
Reading WENO parameters from weno.inp.
Initializing physics. Model = "navierstokes3d"
Reading physical model inputs from file "physics.inp".
Writing solution file iblank_00000.bin.
Setting up time integration.
Setting up libROM interface.
libROM inputs and parameters:
reduced model dimensionality: 32
sampling frequency: 1
mode: train
type: DMD
save to file: true
local vector size: 61440
libROM DMD inputs:
number of samples per window: 100
directory name for DMD onjects: DMD
Solving in time (from 0 to 400 iterations)
Writing immersed body surface data file surface00000.dat.
Writing solution file op_00000.bin.
DMDROMObject::takeSample() - creating new DMD object, t=0.000000 (total: 1).
iter= 10 t=5.000E-02 CFL=5.878E-01 norm=7.4787E-01 wctime: 1.7E+00 (s)
iter= 20 t=1.000E-01 CFL=5.862E-01 norm=6.8901E-01 wctime: 1.9E+00 (s)
iter= 30 t=1.500E-01 CFL=5.859E-01 norm=7.4701E-01 wctime: 2.1E+00 (s)
iter= 40 t=2.000E-01 CFL=5.857E-01 norm=6.8995E-01 wctime: 2.0E+00 (s)
Writing immersed body surface data file surface00001.dat.
Writing solution file op_00001.bin.
iter= 50 t=2.500E-01 CFL=5.856E-01 norm=7.5315E-01 wctime: 2.1E+00 (s)
iter= 60 t=3.000E-01 CFL=1.052E+00 norm=6.7415E-01 wctime: 2.1E+00 (s)
iter= 70 t=3.500E-01 CFL=9.818E-01 norm=7.0618E-01 wctime: 2.1E+00 (s)
iter= 80 t=4.000E-01 CFL=9.227E-01 norm=6.5493E-01 wctime: 2.1E+00 (s)
Writing immersed body surface data file surface00002.dat.
Writing solution file op_00002.bin.
iter= 90 t=4.500E-01 CFL=1.003E+00 norm=7.0464E-01 wctime: 2.1E+00 (s)
iter= 100 t=5.000E-01 CFL=1.150E+00 norm=6.8393E-01 wctime: 2.1E+00 (s)
DMDROMObject::takeSample() - creating new DMD object, t=0.500000 (total: 2).
iter= 110 t=5.500E-01 CFL=1.235E+00 norm=7.1132E-01 wctime: 2.2E+00 (s)
iter= 120 t=6.000E-01 CFL=1.324E+00 norm=6.9216E-01 wctime: 2.5E+00 (s)
Writing immersed body surface data file surface00003.dat.
Writing solution file op_00003.bin.
iter= 130 t=6.500E-01 CFL=1.358E+00 norm=7.0342E-01 wctime: 2.5E+00 (s)
iter= 140 t=7.000E-01 CFL=1.397E+00 norm=7.0115E-01 wctime: 2.4E+00 (s)
iter= 150 t=7.500E-01 CFL=1.449E+00 norm=6.8685E-01 wctime: 2.1E+00 (s)
iter= 160 t=8.000E-01 CFL=1.460E+00 norm=7.0893E-01 wctime: 2.1E+00 (s)
Writing immersed body surface data file surface00004.dat.
Writing solution file op_00004.bin.
iter= 170 t=8.500E-01 CFL=1.463E+00 norm=6.8381E-01 wctime: 2.1E+00 (s)
iter= 180 t=9.000E-01 CFL=1.464E+00 norm=7.0983E-01 wctime: 2.1E+00 (s)
iter= 190 t=9.500E-01 CFL=1.468E+00 norm=6.7942E-01 wctime: 2.1E+00 (s)
iter= 200 t=1.000E+00 CFL=1.480E+00 norm=7.0922E-01 wctime: 2.2E+00 (s)
Writing immersed body surface data file surface00005.dat.
Writing solution file op_00005.bin.
DMDROMObject::takeSample() - creating new DMD object, t=1.000000 (total: 3).
iter= 210 t=1.050E+00 CFL=1.492E+00 norm=6.7655E-01 wctime: 2.2E+00 (s)
iter= 220 t=1.100E+00 CFL=1.497E+00 norm=7.0630E-01 wctime: 2.2E+00 (s)
iter= 230 t=1.150E+00 CFL=1.506E+00 norm=6.7587E-01 wctime: 2.2E+00 (s)
iter= 240 t=1.200E+00 CFL=1.510E+00 norm=7.0913E-01 wctime: 2.2E+00 (s)
Writing immersed body surface data file surface00006.dat.
Writing solution file op_00006.bin.
iter= 250 t=1.250E+00 CFL=1.513E+00 norm=6.7654E-01 wctime: 2.3E+00 (s)
iter= 260 t=1.300E+00 CFL=1.509E+00 norm=7.0671E-01 wctime: 2.6E+00 (s)
iter= 270 t=1.350E+00 CFL=1.505E+00 norm=6.8107E-01 wctime: 2.7E+00 (s)
iter= 280 t=1.400E+00 CFL=1.503E+00 norm=7.0402E-01 wctime: 2.7E+00 (s)
Writing immersed body surface data file surface00007.dat.
Writing solution file op_00007.bin.
iter= 290 t=1.450E+00 CFL=1.501E+00 norm=6.9051E-01 wctime: 2.2E+00 (s)
iter= 300 t=1.500E+00 CFL=1.496E+00 norm=7.0816E-01 wctime: 2.2E+00 (s)
DMDROMObject::takeSample() - creating new DMD object, t=1.500000 (total: 4).
iter= 310 t=1.550E+00 CFL=1.492E+00 norm=6.9908E-01 wctime: 2.2E+00 (s)
iter= 320 t=1.600E+00 CFL=1.487E+00 norm=7.1111E-01 wctime: 2.3E+00 (s)
Writing immersed body surface data file surface00008.dat.
Writing solution file op_00008.bin.
iter= 330 t=1.650E+00 CFL=1.485E+00 norm=7.0298E-01 wctime: 2.3E+00 (s)
iter= 340 t=1.700E+00 CFL=1.481E+00 norm=7.0506E-01 wctime: 2.3E+00 (s)
iter= 350 t=1.750E+00 CFL=1.478E+00 norm=7.0788E-01 wctime: 2.3E+00 (s)
iter= 360 t=1.800E+00 CFL=1.476E+00 norm=7.0432E-01 wctime: 2.3E+00 (s)
Writing immersed body surface data file surface00009.dat.
Writing solution file op_00009.bin.
iter= 370 t=1.850E+00 CFL=1.472E+00 norm=7.1498E-01 wctime: 2.3E+00 (s)
iter= 380 t=1.900E+00 CFL=1.469E+00 norm=7.0543E-01 wctime: 2.2E+00 (s)
iter= 390 t=1.950E+00 CFL=1.471E+00 norm=7.2006E-01 wctime: 2.1E+00 (s)
iter= 400 t=2.000E+00 CFL=1.469E+00 norm=7.0183E-01 wctime: 1.8E+00 (s)
Completed time integration (Final time: 2.000000), total wctime: 815.048244 (seconds).
Writing immersed body surface data file surface00010.dat.
Writing solution file op_00010.bin.
libROM: Training ROM.
DMDROMObject::train() - training DMD object 0 with 101 samples.
Using 32 basis vectors out of 100.
DMDROMObject::train() - training DMD object 1 with 101 samples.
Using 32 basis vectors out of 100.
DMDROMObject::train() - training DMD object 2 with 101 samples.
Using 32 basis vectors out of 100.
DMDROMObject::train() - training DMD object 3 with 100 samples.
Using 32 basis vectors out of 99.
libROM: total training wallclock time: 155.407039 (seconds).
libROM: Predicting solution at time 0.0000e+00 using ROM.
libROM: wallclock time: 2.914239 (seconds).
Writing solution file op_rom_00000.bin.
libROM: Predicting solution at time 2.0000e-01 using ROM.
libROM: wallclock time: 2.913659 (seconds).
Writing solution file op_rom_00001.bin.
libROM: Predicting solution at time 4.0000e-01 using ROM.
libROM: wallclock time: 2.914199 (seconds).
Writing solution file op_rom_00002.bin.
libROM: Predicting solution at time 6.0000e-01 using ROM.
libROM: wallclock time: 2.979725 (seconds).
Writing solution file op_rom_00003.bin.
libROM: Predicting solution at time 8.0000e-01 using ROM.
libROM: wallclock time: 3.099582 (seconds).
Writing solution file op_rom_00004.bin.
libROM: Predicting solution at time 1.0000e+00 using ROM.
libROM: wallclock time: 2.960964 (seconds).
Writing solution file op_rom_00005.bin.
libROM: Predicting solution at time 1.2000e+00 using ROM.
libROM: wallclock time: 2.940771 (seconds).
Writing solution file op_rom_00006.bin.
libROM: Predicting solution at time 1.4000e+00 using ROM.
libROM: wallclock time: 3.098174 (seconds).
Writing solution file op_rom_00007.bin.
libROM: Predicting solution at time 1.6000e+00 using ROM.
libROM: wallclock time: 2.902298 (seconds).
Writing solution file op_rom_00008.bin.
libROM: Predicting solution at time 1.8000e+00 using ROM.
libROM: wallclock time: 2.954182 (seconds).
Writing solution file op_rom_00009.bin.
libROM: Predicting solution at time 2.0000e+00 using ROM.
libROM: wallclock time: 2.914002 (seconds).
libROM: Calculating diff between PDE and ROM solutions.
Writing solution file op_rom_00010.bin.
libROM: total prediction/query wallclock time: 32.591795 (seconds).
libROMInterface::saveROM() - saving ROM objects.
Saving DMD object with filename root DMD/dmdobj_0000.
Saving DMD object with filename root DMD/dmdobj_0001.
Saving DMD object with filename root DMD/dmdobj_0002.
Saving DMD object with filename root DMD/dmdobj_0003.
Norms of the diff between ROM and PDE solutions for domain 0:
L1 Norm : 6.4931422313752370E-03
L2 Norm : 1.9356397835829150E-02
Linfinity Norm : 9.1818439088445988E-02
Solver runtime (in seconds): 1.0099263970000000E+03
Total runtime (in seconds): 1.0104446960000000E+03
Deallocating arrays.
Finished.