Efficient Transonic Wing Flow Simulation and Shock Capture Using Euler‑PINN

Problem Background and Application Scenario

Transonic compressible flow over an airfoil is a classic benchmark where local Mach numbers exceed unity, generating shock waves that cause abrupt jumps in density, pressure, and velocity, and make wall pressure distribution and lift highly sensitive to shock position. Structured O‑type grids are advantageous because they conform to the geometry, allow natural far‑field conditions, and enable local refinement near shocks.

Governing Equations and Boundary Conditions

The steady two‑dimensional Euler equations are expressed in primitive variables (density rho, velocity components u, v, and pressure p). Strong‑form residuals for continuity, momentum (x and y), and energy are constructed and weighted by cell volume in the computational domain. Boundary conditions include far‑field free‑stream values, no‑penetration wall (zero normal velocity), and periodic matching on the left/right sides of the O‑grid.

Coordinate Transformation and Jacobian Chain Rule

The physical domain is mapped to a rectangular computational domain. Grid indices are normalized (zero mean, unit variance) and used as network inputs. Jacobian matrices are pre‑computed on the grid to convert derivatives from computational to physical coordinates, enabling efficient evaluation of the PDE residuals via JAX automatic differentiation.

Network Architecture and Tensor Flow

A weight‑normalized multilayer perceptron with 5 hidden layers (128 neurons each) takes the 2‑D normalized grid indices and outputs the four primitive variables. The tensor flow from grid indices to PDE residuals is summarized in a table (omitted for brevity). The network is trained on 20 L‑BFGS epochs, each using a random mini‑batch of 16 080 points (≈80 % of the 20 100 grid nodes).

Core Code: PDE Residual Computation

# pred: (batch, 4) → u, v, p, rho
# Dxy: (batch, 2, 4) → physical gradients [∂/∂x, ∂/∂y] × [u,v,p,ρ]

u_x, v_x, p_x, rho_x = Dxy[:,0,0], Dxy[:,0,1], Dxy[:,0,2], Dxy[:,0,3]
u_y, v_y, p_y, rho_y = Dxy[:,1,0], Dxy[:,1,1], Dxy[:,1,2], Dxy[:,1,3]

a2 = GAMMA * p / rho

r1 = rho*u_x + u*rho_x + rho*v_y + v*rho_y          # continuity
r2 = u*u_x + v*u_y + p_x/rho                       # x‑momentum
r3 = u*v_x + v*v_y + p_y/rho                       # y‑momentum
r4 = rho*a2*(u_x+v_y) + u*p_x + v*p_y               # energy

w = Jsq[batch_idx]                                 # volume weight
loss_pde = ((r1**2 + r2**2 + r3**2 + r4**2) * w).sum() / w.sum()

Loss Function and Optimization Strategy

The total loss combines weighted PDE residuals and three boundary‑loss terms (far‑field, wall, periodic) with comparable magnitudes at the start of training. L‑BFGS (memory = 100) performs up to 1 000 inner iterations per epoch; the training wall‑time is 55.3 s. Safety mechanisms automatically reset the random seed if the loss becomes non‑finite and roll back parameters when a loss spike exceeds four times the previous epoch.

Reference Solution and Accuracy Assessment

A high‑resolution finite‑volume method (MUSCL reconstruction, minmod limiter, HLLC Riemann solver, third‑order TVD Runge‑Kutta, CFL = 0.9) is solved on the same O‑grid to provide a full‑field reference ( fvm_reference_field_200x100.npz) and a denser wall‑pressure reference ( Airfoil4_800_400_Cp.dat). Relative errors are:

Pressure: 0.65 %

Velocity magnitude: 0.97 %

Density: 0.52 %

Normal velocity component near the shock: 10.26 %

Wall‑pressure coefficient: 0.136 % (relative)

Lift coefficient: 15.5 % (relative)

Errors concentrate in the shock region; the smooth PINN struggles to represent the discontinuity, leading to larger relative errors for the normal‑velocity component.

Efficiency vs. Accuracy Trade‑offs

Compared with the FVM reference, Euler‑PINN‑JAX requires only 55 s of wall‑time versus a full FVM solve, but the strong‑form smoothness constraint limits shock‑capture accuracy. Potential improvements include weak‑form residuals, explicit shock‑location parameterization, or adaptive weighting of points in high‑gradient regions.

Common Questions and Solutions

Q1: Low PDE loss does not guarantee low overall error because the loss is evaluated on a random subset and shock regions occupy a tiny area.

Q2: The normal‑velocity error is amplified by the weak discontinuity across the airfoil surfaces and the steep gradient of the shock.

Q3: Sudden loss spikes are mitigated by an automatic rollback mechanism.

Q4: Increasing point density near the shock helps but cannot fully overcome the smoothness limitation of strong‑form PINNs.

Q5: JAX provides just‑in‑time compilation ( jit) and efficient Jacobian‑vector products ( jvp) that accelerate L‑BFGS on GPUs.

Conclusion and Outlook

The tutorial demonstrates a complete workflow for solving transonic airfoil flow with Euler‑PINN in JAX, achieving sub‑minute training times and reasonable predictions in smooth regions, while highlighting intrinsic limitations in shock capture. Future work will explore weak‑form formulations, learnable shock interfaces, and adaptive sampling to improve accuracy for discontinuous solutions.