See 1D Euler Equations - Sod Shock Tube to familiarize yourself with this case. This example uses a DMD object that has already been trained (see 1D Euler Equations - Sod Shock Tube (Time Windowed DMD)).

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

Governing equations: 1D Euler equations (euler1d.h)

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.

References:

G.A. Sod, "A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws," J. Comput. Phys., 27, 1 (1978).
C. B. Laney, "Computational Gasdynamics", Cambridge University Press, 1998.

Domain: \(0 \le x \le 1.0\), "extrapolate" (_EXTRAPOLATE_) boundary conditions

Initial Solution:

\( 0 \le x < 0.5\): \(\rho = 1, u = 0, p = 1\)
\( 0.5 \le x \le 1\): \(\rho = 0.125, u = 0, p = 0.1\)

Reduced Order Modeling:

Mode: predict (libROMInterface::m_mode, _ROM_MODE_PREDICT_)
Component Mode: monolithic (libROMInterface::m_comp_mode, _ROM_COMP_MODE_MONOLITHIC_)
Type: Dynamic Mode Decomposition (DMD) with time windowing (libROMInterface::m_rom_type)

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 1D Euler Equations - Sod Shock Tube (Time Windowed DMD).

solver.inp:

begin
  ndims             1
  nvars             3
  size              201
  iproc             10
  ghost             3
  n_iter            20
  dt                0.01
  screen_op_iter    1
  file_op_iter      1
  op_file_format    binary
  op_overwrite      no
  model             euler1d
end

boundary.inp

2
extrapolate      0     1     0     0
extrapolate      0    -1     0     0

physics.inp

begin
  gamma       1.4
  upwinding   rf-char
end

To generate initial.inp, compile and run the following code in the run directory:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
int main(){
  
    int NI=101,ndims=1;
  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) {
    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 { 
      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-1));
    double *x, *rho,*rhou,*e;
    x    = (double*) calloc (NI, sizeof(double));
    rho  = (double*) calloc (NI, sizeof(double));
    rhou = (double*) calloc (NI, sizeof(double));
    e    = (double*) calloc (NI, sizeof(double));
    for (i = 0; i < NI; i++){
        x[i] = i*dx;
    double RHO,U,P;
    if (x[i] < 0.5) {
      RHO = 1.0;
      U   = 0.0;
      P   = 1.0;
    } else {
      RHO = 0.125;
      U   = 0;
      P   = 0.1;
    }
    rho[i]  = RHO;
    rhou[i] = RHO*U;
    e[i]    = P/0.4 + 0.5*RHO*U*U;
    }
  if (!strcmp(ip_file_type,"ascii")) {
    FILE *out;
      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 ",rho[i]);                     
    fprintf(out,"\n");
      for (i = 0; i < NI; i++)  fprintf(out,"%lf ",rhou[i]);                        
    fprintf(out,"\n");
      for (i = 0; i < NI; i++)  fprintf(out,"%lf ",e[i]);                       
    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(rho);
    free(rhou);
    free(e);
    return(0);
}

Output:

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

21 output files op_00000.bin, op_00001.bin, ... op_00020.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=0.2\). 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_1DBinary.py) can be used to generate plots from these binary files. Alternatively, HyPar::op_file_format can be set to text, and GNUPlot or something similar can be used to plot the resulting text files.

The animation shows the evolution of the solution (density).

Expected screen output:

HyPar - Parallel (MPI) version with 10 processes
Compiled with PETSc time integration.
Allocated simulation object(s).
Reading solver inputs from file "solver.inp".
  No. of dimensions                          : 1
  No. of variables                           : 3
  Domain size                                : 201 
  Processes along each dimension             : 10 
  Exact solution domain size                 : 201 
  No. of ghosts pts                          : 3
  No. of iter.                               : 20
  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-02
  Check for conservation                     : no
  Screen output iterations                   : 1
  File output iterations                     : 1
  Initial solution file type                 : ascii
  Initial solution read mode                 : serial
  Solution file write mode                   : serial
  Solution file format                       : binary
  Overwrite solution file                    : no
  Physical model                             : euler1d
Partitioning domain and allocating data arrays.
Reading array from ASCII file initial.inp (Serial mode).
Volume integral of the initial solution:
 0:  5.6312499999999999E-01
 1:  0.0000000000000000E+00
 2:  1.3762500000000000E+00
Reading boundary conditions from boundary.inp.
  Boundary                    extrapolate:  Along dimension  0 and face +1
  Boundary                    extrapolate:  Along dimension  0 and face -1
2 boundary condition(s) read.
Initializing solvers.
Initializing physics. Model = "euler1d"
Reading physical model inputs from file "physics.inp".
Setting up libROM interface.
libROM inputs and parameters:
  reduced model dimensionality:  17860400
  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,2.00e-02].
  Loading DMD object (DMD/dmdobj_0001), time window=[2.00e-02,4.00e-02].
  Loading DMD object (DMD/dmdobj_0002), time window=[4.00e-02,6.00e-02].
  Loading DMD object (DMD/dmdobj_0003), time window=[6.00e-02,8.00e-02].
  Loading DMD object (DMD/dmdobj_0004), time window=[8.00e-02,1.00e-01].
  Loading DMD object (DMD/dmdobj_0005), time window=[1.00e-01,1.20e-01].
  Loading DMD object (DMD/dmdobj_0006), time window=[1.20e-01,1.40e-01].
  Loading DMD object (DMD/dmdobj_0007), time window=[1.40e-01,1.60e-01].
  Loading DMD object (DMD/dmdobj_0008), time window=[1.60e-01,1.80e-01].
  Loading DMD object (DMD/dmdobj_0009), time window=[1.80e-01,-1.00e+00].
libROM: Predicted solution at time 0.0000e+00 using ROM, wallclock time: 0.000227.
Writing solution file op_00000.bin.
libROM: Predicted solution at time 1.0000e-02 using ROM, wallclock time: 0.089848.
Writing solution file op_00001.bin.
libROM: Predicted solution at time 2.0000e-02 using ROM, wallclock time: 0.078646.
Writing solution file op_00002.bin.
libROM: Predicted solution at time 3.0000e-02 using ROM, wallclock time: 0.109507.
Writing solution file op_00003.bin.
libROM: Predicted solution at time 4.0000e-02 using ROM, wallclock time: 0.103133.
Writing solution file op_00004.bin.
libROM: Predicted solution at time 5.0000e-02 using ROM, wallclock time: 0.170235.
Writing solution file op_00005.bin.
libROM: Predicted solution at time 6.0000e-02 using ROM, wallclock time: 0.084616.
Writing solution file op_00006.bin.
libROM: Predicted solution at time 7.0000e-02 using ROM, wallclock time: 0.060489.
Writing solution file op_00007.bin.
libROM: Predicted solution at time 8.0000e-02 using ROM, wallclock time: 0.104553.
Writing solution file op_00008.bin.
libROM: Predicted solution at time 9.0000e-02 using ROM, wallclock time: 0.034927.
Writing solution file op_00009.bin.
libROM: Predicted solution at time 1.0000e-01 using ROM, wallclock time: 0.084440.
Writing solution file op_00010.bin.
libROM: Predicted solution at time 1.1000e-01 using ROM, wallclock time: 0.087856.
Writing solution file op_00011.bin.
libROM: Predicted solution at time 1.2000e-01 using ROM, wallclock time: 0.033666.
Writing solution file op_00012.bin.
libROM: Predicted solution at time 1.3000e-01 using ROM, wallclock time: 0.098264.
Writing solution file op_00013.bin.
libROM: Predicted solution at time 1.4000e-01 using ROM, wallclock time: 0.069375.
Writing solution file op_00014.bin.
libROM: Predicted solution at time 1.5000e-01 using ROM, wallclock time: 0.096944.
Writing solution file op_00015.bin.
libROM: Predicted solution at time 1.6000e-01 using ROM, wallclock time: 0.132303.
Writing solution file op_00016.bin.
libROM: Predicted solution at time 1.7000e-01 using ROM, wallclock time: 0.068642.
Writing solution file op_00017.bin.
libROM: Predicted solution at time 1.8000e-01 using ROM, wallclock time: 0.085706.
Writing solution file op_00018.bin.
libROM: Predicted solution at time 1.9000e-01 using ROM, wallclock time: 0.064509.
Writing solution file op_00019.bin.
libROM: Predicted solution at time 2.0000e-01 using ROM, wallclock time: 0.075884.
Writing solution file op_00020.bin.
libROM: total prediction/query wallclock time: 1.733770 (seconds).
Solver runtime (in seconds): 2.9325660000000000E+00
Total  runtime (in seconds): 3.0228090000000001E+00
Deallocating arrays.
Finished.