Having problems computing PDE residuals for The_Well datasets (e.g., shear_flow)

[LeslisXu]

Hi all — I'm computing PDE residuals for The_Well datasets (e.g. turbulent_radiative_layer_2D and shear_flow) using finite differences, but the residuals are much larger than I expect. The data are generated by simulations, so I expected residuals close to zero.

What I do

1. I compute derivatives with a central difference scheme (second-order in the interior, one-sided at boundaries). Here is my derivative function:

    import torch
   
    def compute_derivative(u, delta, order=1, axis=1):
        """
        Compute numerical derivatives along a specified axis.
   
        Args:
            u: Input tensor with shape [B, T, X, Y] (or [B, T, X, Y, C] for vector fields).
            delta: grid spacing along the chosen axis.
            order: 1 for first derivative, 2 for second derivative.
            axis: 1=time, 2=x, 3=y.
        """
        derivative = torch.zeros_like(u)
   
        if order == 1:
            # second-order central differences in interior, one-sided at boundaries
            if axis == 1:  # time
                derivative[:, 1:-1, :, :] = (u[:, 2:, :, :] - u[:, :-2, :, :]) / (2.0 * delta)
                derivative[:, 0:1, :, :]  = (-3.0*u[:, 0:1, :, :] + 4.0*u[:, 1:2, :, :] - u[:, 2:3, :, :]) / (2.0 * delta)
                derivative[:, -1:, :, :]  = (3.0*u[:, -1:, :, :] - 4.0*u[:, -2:-1, :, :] + u[:, -3:-2, :, :]) / (2.0 * delta)
            elif axis == 2:  # x
                derivative[:, :, 1:-1, :] = (u[:, :, 2:, :] - u[:, :, :-2, :]) / (2.0 * delta)
                derivative[:, :, 0:1, :]  = (-3.0*u[:, :, 0:1, :] + 4.0*u[:, :, 1:2, :] - u[:, :, 2:3, :]) / (2.0 * delta)
                derivative[:, :, -1:, :]  = (3.0*u[:, :, -1:, :] - 4.0*u[:, :, -2:-1, :] + u[:, :, -3:-2, :]) / (2.0 * delta)
            elif axis == 3:  # y
                derivative[:, :, :, 1:-1] = (u[:, :, :, 2:] - u[:, :, :, :-2]) / (2.0 * delta)
                derivative[:, :, :, 0:1]  = (-3.0*u[:, :, :, 0:1] + 4.0*u[:, :, :, 1:2] - u[:, :, :, 2:3]) / (2.0 * delta)
                derivative[:, :, :, -1:]  = (3.0*u[:, :, :, -1:] - 4.0*u[:, :, :, -2:-1] + u[:, :, :, -3:-2]) / (2.0 * delta)
   
        elif order == 2:
            # second-order second derivative
            if axis == 1:
                derivative[:, 1:-1, :, :] = (u[:, 2:, :, :] - 2*u[:, 1:-1, :, :] + u[:, :-2, :, :]) / (delta**2)
                derivative[:, 0:1, :, :]  = (u[:, 2:3, :, :] - 2*u[:, 1:2, :, :] + u[:, 0:1, :, :]) / (delta**2)
                derivative[:, -1:, :, :]  = (u[:, -1:, :, :] - 2*u[:, -2:-1, :, :] + u[:, -3:-2, :, :]) / (delta**2)
            elif axis == 2:
                derivative[:, :, 1:-1, :] = (u[:, :, 2:, :] - 2*u[:, :, 1:-1, :] + u[:, :, :-2, :]) / (delta**2)
                derivative[:, :, 0:1, :]  = (u[:, :, 2:3, :] - 2*u[:, :, 1:2, :] + u[:, :, 0:1, :]) / (delta**2)
                derivative[:, :, -1:, :]  = (u[:, :, -1:, :] - 2*u[:, :, -2:-1, :] + u[:, :, -3:-2, :]) / (delta**2)
            elif axis == 3:
                derivative[:, :, :, 1:-1] = (u[:, :, :, 2:] - 2*u[:, :, :, 1:-1] + u[:, :, :, :-2]) / (delta**2)
                derivative[:, :, :, 0:1]  = (u[:, :, :, 2:3] - 2*u[:, :, :, 1:2] + u[:, :, :, 0:1]) / (delta**2)
                derivative[:, :, :, -1:]  = (u[:, :, :, -1:] - 2*u[:, :, :, -2:-1] + u[:, :, :, -3:-2]) / (delta**2)
   
        return derivative
   

2. From the HDF5 file I load fields with shapes:

    pressure: (4, 200, 256, 512)
    tracer:   (4, 200, 256, 512)
    velocity: (4, 200, 256, 512, 2)
   

   I call compute_physics_loss(...) with velocity[...,0] and velocity[...,1] as velocity_x and velocity_y.

3. The PDEs I use:

   Momentum (component form)

   $$ \partial_t \mathbf{u} + \nabla p - \nu \Delta \mathbf{u} = -\mathbf{u}\cdot\nabla\mathbf{u}. $$

   Tracer

   $$ \partial_t s - D\Delta s = -\mathbf{u}\cdot\nabla s, $$

   where ν = 1/Reynolds and D = ν/Schmidt.

   and the PDE terms are calculated using

    # time derivatives (one-order)
    du_x_dt = compute_derivative(velocity_x, dt, order=1, axis=1)
    du_y_dt = compute_derivative(velocity_y, dt, order=1, axis=1)
    ds_dt = compute_derivative(tracer, dt, order=1, axis=1)
    
    # spatial derivatives
    # pressure gradients (one-order)
    dp_dx = compute_derivative(pressure, dx, order=1, axis=2)
    dp_dy = compute_derivative(pressure, dy, order=1, axis=3)
   
    # velocity Laplacian (second-order)
    du_x_dx2 = compute_derivative(velocity_x, dx, order=2, axis=2)
    du_x_dy2 = compute_derivative(velocity_x, dy, order=2, axis=3)
    laplacian_u_x = du_x_dx2 + du_x_dy2
    
    du_y_dx2 = compute_derivative(velocity_y, dx, order=2, axis=2)
    du_y_dy2 = compute_derivative(velocity_y, dy, order=2, axis=3)
    laplacian_u_y = du_y_dx2 + du_y_dy2
    
    # tracer Laplacian (second-order)
    ds_dx2 = compute_derivative(tracer, dx, order=2, axis=2)
    ds_dy2 = compute_derivative(tracer, dy, order=2, axis=3)
    laplacian_s = ds_dx2 + ds_dy2
    
    # convection term (u·∇u)
    du_x_dx = compute_derivative(velocity_x, dx, order=1, axis=2)
    du_x_dy = compute_derivative(velocity_x, dy, order=1, axis=3)
    advection_u_x = velocity_x * du_x_dx + velocity_y * du_x_dy
    
    du_y_dx = compute_derivative(velocity_y, dx, order=1, axis=2)
    du_y_dy = compute_derivative(velocity_y, dy, order=1, axis=3)
    advection_u_y = velocity_x * du_y_dx + velocity_y * du_y_dy
   
    # convection term (u·∇s)
    ds_dx = compute_derivative(tracer, dx, order=1, axis=2)
    ds_dy = compute_derivative(tracer, dy, order=1, axis=3)
    advection_s = velocity_x * ds_dx + velocity_y * ds_dy
   

4. I compute residual tensors (shape [B, T, X, Y]) as:

    momentum_x_residual = du_x_dt + dp_dx - nu * laplacian_u_x + advection_u_x
    momentum_y_residual = du_y_dt + dp_dy - nu * laplacian_u_y + advection_u_y
    tracer_residual     = ds_dt - D * laplacian_s + advection_s
   

5. I then compute the mean squared residuals via:

    res_mom_x = ((momentum_x_residual)**2).mean().item()
    res_mom_y = ((momentum_y_residual)**2).mean().item()
    res_tracer = ((tracer_residual)**2).mean().item()
   

   (Yes, that expression is the mean squared residual / MSE.)

Observed output

The momentum equation residual norms: 0.18399687111377716, 0.1050143912434578
The tracer equation residual norm: 2.9967291355133057

I expected values near zero but get values ~0.1–3.0.

What I checked so far

• data shapes and indexing
• derivative implementation (central diff interior + one-sided boundary)

Questions / help requested

• Am I missing any common pitfalls that could cause such large residuals?
• Could the dataset use different conventions (e.g. staggered grid, pressure vs velocity staggering, units/scaling, periodic BCs, ghost cells, or coordinate transforms) that I should account for?
• Are there mistakes in my finite-difference discretization or sign conventions for the PDE terms?
• Are there recommended sanity checks I should run to localize the issue?

Relative Information

dt=0.09999999999999999, dx=0.003921568859368563, dy=0.001956947147846222, Reynolds=500000.0, Schmidt=0.20000000298023224

Thanks in advance for any pointers!!!

[mikemccabe210]

Hi @LeslisXu,

I haven't checked your specific implementation, but I wouldn't expect finite-differences to be a good approximation for most data in the Well or for non-trivial simulation data in general unless the recorded time intervals are extremely small. In the case of most Well data:

1. Time sampling in most datasets is significantly sparser than the base integrator. Generally the underlying solver is using a CFL-based time step criteria for most of the CFD-oriented datasets. TRL2D for instance has a typical simulation time step of O(10^-2) while it's recorded at O(0). I don't know off-hand what the simulation step for shear-flow was, but the timestep is computed using a second order IMEX method and likely follows a similar ratio between integrator steps and recorded steps.
2. Very few datasets were generated with FD spatial discretizations. For the two examples, TRL uses athena++ which I believe is fundamentally finite-volume code with AMR while shear flow is using a pseudospectral discretization.

So the residual is computing the difference in endpoints between a solver utilizing low order FD in space and explicit Euler in time with possibly nonlinear higher order in time and space methods using significantly smaller steps over a fixed time horizon which generally won't give you similar answers unless both discretizations are stable and the system is converging to a fixed point.

Most solver details can be found on the datasets page, though different contributors providing different levels of detail there.

Best,
Mike

[LeslisXu]

Hi @mikemccabe210 ,
Thank you for the detailed explanation! This is really helpful.
I understand now—the sparse temporal sampling and the differences between the original high-order methods (finite volume, pseudospectral) and simple finite differences make FD approximations unsuitable for computing accurate residuals on this data.
I'll check the datasets page for the specific solver details and reconsider my approach accordingly.
Thank you again for clarifying this!
Best,
Leslie