- Я вычисляю производные с помощью центральной разностной схемы (второй порядок внутри, односторонний на границах). Вот моя производная функция:
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 загружаю поля с фигурами:
pressure: (4, 200, 256, 512)
tracer: (4, 200, 256, 512)
velocity: (4, 200, 256, 512, 2)
Я вызываю Compute_Physics_loss(...) со скоростью[...,0] и Velocity[...,1] как Velocity_x и Velocity_y. - Используемые мной PDE:
Импульс (форма компонента)
$$ \partial_t \mathbf{u} + \nabla p - \nu \Delta \mathbf{u} = -\mathbf{u}\cdot\nabla\mathbf{u}. $$
Трейсер
$$ \partial_t s - D\Delta s = -\mathbf{u}\cdot\nabla s, $$
где ν = 1/Рейнольдс и D = ν/Шмидт.
и члены 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 - Затем я вычисляю среднеквадратические остатки с помощью:
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()
(Да, это выражение представляет собой среднеквадратический остаток / MSE.)
- Нормы остатка уравнения импульса: 0,18399687111377716, 0,1050143912434578
- Норма остатка уравнения трассировки: 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
