Интегрирование векторного поля в трехмерном пространстве (как обратное градиенту) ⇐ Python

[Anonymous]

Я хотел бы написать функцию (на Python), которая интегрирует векторное поле в трехмерную сетку,
скажем, обратную градиенту.
Я видел несколько примеров и логика проста, но я получаю неточные численные результаты, когда проверяю их, снова реконструируя градиент с помощью numpy.gradient.
Как бы мне это изменить? С самой процедурой все в порядке?
Вот пример кода, который я пробовал:

Код: Выделить всё

#!/usr/bin/env python
# assume hypothetical smooth 3d vector field
# grid 51x51x51
# vector component smooth changing -2~2

import numpy as np
from scipy.ndimage import gaussian_filter

## This block: generate a random smooth data of vector field in 3D space
Nx, Ny, Nz = 51, 51, 51
dx, dy, dz = 0.1, 0.1, 0.1
min_val, max_val = -2, 2
vector_field_3d = np.random.uniform(min_val, max_val, (Nx, Ny, Nz, 3))
for i in range(3):
vector_field_3d[...,i] = gaussian_filter(vector_field_3d[...,i], sigma=3)
max_component = np.max(np.abs(vector_field_3d))
if max_component >  0:
vector_field_3d = vector_field_3d / max_component * (max_val - min_val) / 2
vector_field_table = vector_field_3d.reshape(-1,3)

np.savetxt('smoothtest_vector_field_table.dat',vector_field_table,fmt='%15.8f')

## Try integration, and validate the result with reconstructing gradient again
def reverse_gradient_integrate_tripleloop(table3d_vector_field,grid_spacing):
dx, dy, dz = grid_spacing
table3d_integ = np.zeros(table3d_vector_field.shape[:3])
#dxintegral = np.cumsum(table3d_vector_field[:,:,:,0], axis=0)*dx
#dyintegral = np.cumsum(table3d_vector_field[:,:,:,1], axis=1)*dy
#dzintegral = np.cumsum(table3d_vector_field[:,:,:,2], axis=2)*dz
nx, ny, nz, nq = table3d_vector_field.shape
for i in range(1,nx):
table3d_integ[i,0,0] = dx*(table3d_vector_field[i-1,0,0,0]) + table3d_integ[i-1,0,0]
for j in range(1,ny):
table3d_integ[i,j,0] = dy*(table3d_vector_field[i,j-1,0,1]) + table3d_integ[i,j-1,0]
for k in range(1,nz):
table3d_integ[i,j,k] = dz*(table3d_vector_field[i,j,k-1,2]) + table3d_integ[i,j,k-1]
return table3d_integ

def reverse_gradient_integrate_diagonal(table3d_vector_field,grid_spacing):
dx, dy, dz = grid_spacing
fx, fy, fz = table3d_vector_field[:,:,:,0], table3d_vector_field[:,:,:,1], table3d_vector_field[:,:,:,2]
nx, ny, nz, nq = table3d_vector_field.shape
table3d_integ = np.zeros(table3d_vector_field.shape[:3])
# diagonal consideration: for every grid, take all 3 f_xdx,f_ydy,f_zdz values adjacent by [i-1,j-1,k-1] cell
# this way uses full information.
# first, complete (1) edges (2) surfaces (2d diagonal from edge values) then (3) inner values
# (1) 3 edges
table3d_integ[:,0,0] = np.cumsum( fx[:,0,0] )*dx
table3d_integ[0,:,0] = np.cumsum( fy[0,:,0] )*dy
table3d_integ[0,0,:] = np.cumsum( fz[0,0,:] )*dz
# (2) surfaces
for i in range(1,nx):
for j in range(1,ny):
table3d_integ[i,j,0] = dx*fx[i-1,j-1,0] + dy*fy[i-1,j-1,0] + table3d_integ[i-1,j-1,0]
for i in range(1,nx):
for k in range(1,nz):
table3d_integ[i,0,k] = dx*fx[i-1,0,k-1] + dz*fz[i-1,0,k-1] + table3d_integ[i-1,0,k-1]
for j in range(1,ny):
for k in range(1,nz):
table3d_integ[0,j,k] = dy*fy[0,j-1,k-1] + dz*fz[0,j-1,k-1] + table3d_integ[0,j-1,k-1]
# (3) full inner elements
for i in range(1,nx):
for j in range(1,ny):
for k in range(1,nz):
table3d_integ[i,j,k] = dx*fx[i-1,j-1,k-1] + dy*fy[i-1,j-1,k-1] + dz*fz[i-1,j-1,k-1] + table3d_integ[i-1,j-1,k-1]
return table3d_integ

def reverse_gradient_integrate_doublediagonal(table3d_vector_field,grid_spacing):
table3d_integ_forward = reverse_gradient_integrate_diagonal(table3d_vector_field,grid_spacing)
table3d_vector_field_rev = table3d_vector_field[::-1,::-1,::-1,:]
table3d_vector_field_revm = -1.0*table3d_vector_field_rev
table3d_integ_backward_rev = reverse_gradient_integrate_diagonal(table3d_vector_field_revm,grid_spacing)
table3d_integ_backward = table3d_integ_backward_rev[::-1,::-1,::-1]
table3d_integ = (table3d_integ_forward + table3d_integ_backward) / 2.0
return table3d_integ

nx, ny, nz = 51, 51, 51
dx1, dy1, dz1 = 0.1, 0.1, 0.1

force_q3_columns = np.loadtxt('smoothtest_vector_field_table.dat')

force_q3 = force_q3_columns.reshape(nx,ny,nz,3)

#table1_energy3d = reverse_gradient_integrate_tripleloop(force_q3,(dx1,dy1,dz1))
table1_energy3d = reverse_gradient_integrate_doublediagonal(force_q3,(dx1,dy1,dz1))
table1_energy3d -= table1_energy3d[-1,-1,-1]
table1_energy3d *= -1.0

## Appendix: gradient calculation again for validation
gf3d_x, gf3d_y, gf3d_z = np.gradient(table1_energy3d,dx1,dy1,dz1,edge_order=1)
gf3d_x *= -1.0
gf3d_y *= -1.0
gf3d_z *= -1.0
recon_force_q3 = np.stack((gf3d_x, gf3d_y, gf3d_z), axis=-1)  # (nx, ny, nz, 3) array
gradient_recon_error = recon_force_q3 - force_q3
fmag_recon_error = np.linalg.norm(gradient_recon_error,axis=-1)
print('Mean latent gradient force reconstruction error check: ',np.mean(fmag_recon_error))

print('reconstructed force q3 grid data shape:', recon_force_q3.shape)
recon_force_table = recon_force_q3.reshape(-1,3)
np.savetxt('reconsmoothtest_vector_field_table.dat',recon_force_table,fmt='%15.8f')

Мое выходное сообщение сообщает следующее:
Проверка ошибки реконструкции средней скрытой градиентной силы: 1,4047385997086408
Что является существенно большой нормой средней ошибки, учитывая величины исходного образца.
Прикрепляя блок кода, я сначала попробовал простейшую прямую интеграцию (

Код: Выделить всё

reverse_gradient_integrate_tripleloop

) и получил большую числовую ошибку,
поэтому переключился на функцию диагонального вычисления интегралов с учетом всех приращений трех осей fxdx,fydy,fzdz из (i-1,j-1,k-1) ячейка.
Затем, учитывая, что np.gradient считает данные как прямой, так и обратной ячейки для расчета наклона, я попытался усреднить два результата прямой и обратной интеграции.
/>

Код: Выделить всё

reverse_gradient_integrate_doublediagonal

несколько уменьшил ошибку, но она все равно большая.
Является ли эта проблема по своей сути неуникальной? или вы бы указали, что функции, которые я использую, недействительны?
Я буду признателен за любой совет. Спасибо!

Подробнее здесь: https://stackoverflow.com/questions/791 ... f-gradient