Что я делаю
- Я вычисляю производные с помощью центральной разностной схемы (второй порядок внутри, односторонний на границах). Вот моя производная функция:
Код: Выделить всё
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 - Из файла HDF5 загружаю поля с фигурами:
Я вызываю Compute_Physics_loss(...) со скоростью[...,0] и Velocity[...,1] как Velocity_x и Velocity_y.
Код: Выделить всё
pressure: (4, 200, 256, 512) tracer: (4, 200, 256, 512) velocity: (4, 200, 256, 512, 2) - Используемые мной PDE:
Импульс (форма компонента)
$$ \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, $$
где ν = 1/Reynolds и D = ν/Schmidt.
и вычисляются члены PDE используяКод: Выделить всё
# 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 - Я вычисляю тензоры остатков (форма [B, T, X, Y]) как:
Код: Выделить всё
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 - Затем я вычисляю среднеквадратические остатки с помощью:
(Да, это выражение представляет собой среднеквадратичный остаток / MSE.)
Код: Выделить всё
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()
Код: Выделить всё
The momentum equation residual norms: 0.18399687111377716, 0.1050143912434578
The tracer equation residual norm: 2.9967291355133057
То, что я проверял на данный момент
- формы данных и индексирование
- производная реализация (центральная внутренняя разница + односторонняя граница)
- Упускаю ли я какие-либо распространенные ошибки, которые могут привести к таким большим остаткам?
- Может ли набор данных использовать другие соглашения (например, шахматную сетку, сдвиг давления в зависимости от скорости, единицы измерения/масштабирование, периодические BC, призрачные ячейки или преобразования координат), которые мне следует учитывать?
- Есть ли ошибки в моей конечно-разностной дискретизации или соглашениях о знаках для терминов PDE?
- Существуют ли рекомендуемые проверки работоспособности, которые мне следует выполнить, чтобы локализовать проблему?
Код: Выделить всё
dt=0.09999999999999999, dx=0.003921568859368563, dy=0.001956947147846222, Reynolds=500000.0, Schmidt=0.20000000298023224
Подробнее здесь: https://stackoverflow.com/questions/797 ... -residuals
