Расчета электрического поля заряженного кольца


Этот скрипт вычисляет формы электрического поля на 2D сетке перпендикулярно равномерно заряженных кольца, используя составляющей по dblquad (руководство, документация).

Я выбрал грубой сетке и указан очень большой относительной погрешности epsrel вот так, что этот пример будет работать около 10 секунд, но я буду использовать его на очень мелкую сетку и с гораздо ниже указанной ошибки.

Я хотел бы знать, если есть способы это может быть значительно ускорена без потери точности, которая должна быть в 1Е-7 микрорайоне.

Я буду избегать патологических точек вручную, как я сделал здесь, избегая плоскости z=0. Я ушла из констант и единиц, чтобы уменьшить беспорядок.

enter image description here

def Excalc(r, th):
    x, y, z = r*np.cos(th), r*np.sin(th), 0.0
    return (x0-x) * ((x0-x)**2 + (y0-y)**2 + (z0-z)**2)**-1.5

def Eycalc(r, th):
    x, y, z = r*np.cos(th), r*np.sin(th), 0.0
    return (y0-y) * ((x0-x)**2 + (y0-y)**2 + (z0-z)**2)**-1.5

def Ezcalc(r, th):
    x, y, z = r*np.cos(th), r*np.sin(th), 0.0
    return (z0-z) * ((x0-x)**2 + (y0-y)**2 + (z0-z)**2)**-1.5


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad

twopi = 2.*np.pi

# annulus of uniform, unit charge density
rmin, rmax   = 0.8, 1.2
thmin, thmax = 0,   twopi

x = np.arange(-2,    2.1, 0.1)
z = np.arange(-1.55, 1.6, 0.1) # avoid z=0 plane for integration.
X, Z = np.meshgrid(x, z)
s = X.shape

eps = 1E-4
y0  = 0

fields, errors = [], []
for x0, z0 in zip(X.flatten(), Z.flatten()):
    # use of lambda in lieu of function calls https://stackoverflow.com/a/49441680/3904031
    Ex, Exerr  = dblquad(Excalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
    Ey, Eyerr  = dblquad(Eycalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
    Ez, Ezerr  = dblquad(Ezcalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
    fields.append([Ex, Ey, Ez])
    errors.append([Exerr, Eyerr, Ezerr])

Evec = np.stack([np.array(thing).reshape(s) for thing in zip(*fields)])
Emag = np.sqrt((Evec**2).sum(axis=0))

Evecerr = np.stack([np.array(thing).reshape(s) for thing in zip(*errors)])
Emagerr = np.sqrt((Evecerr**2).sum(axis=0))

if True:
    names = 'Ex', 'Ey', 'Ez', 'Emag'
    plt.figure()
    for i, (E, name) in enumerate(zip([e for e in Evec] + [Emag], names)):
        smalls = np.abs(E) < 1E-08
        E[smalls] = 0.     # fudge to keep Ey from showing 1E-14 fields
        plt.subplot(2, 2, i+1)
        plt.imshow(E, origin='lower')
        plt.colorbar()
        plt.title(name, fontsize=16)
    plt.show()

if True:
    names = 'Ex rel error', 'Ey rel error', 'Ez rel error', 'Emag rel error'
    plt.figure()
    for i, (err, name) in enumerate(zip([er for er in Evecerr] + [Emagerr], names)):
        plt.subplot(2, 2, i+1)
        plt.imshow(err/E, origin='lower')
        plt.colorbar()
        plt.title(name, fontsize=16)
    plt.show()


343
8
задан 23 марта 2018 в 03:03 Источник Поделиться
Комментарии
2 ответа

Я думаю, что главная трудоемким фактором здесь является расчет dbquadсказал, что есть небольшие улучшения, которые вы можете сделать

разбить на функции

просмотра локальных переменных в функции быстрее , чем глобальный поиск, поэтому putthing каждой части в функции может ускорить этот процесс. Если вы сделаете это, будьте осторожны, что ваш E_calc функции использования глобального состояния (x0, y0 и z0), поэтому они должны быть определены в функции объема, поэтому они могут использовать локальные переменные

и NumPy

Вы можете использовать numpy немного больше, чтобы собрать окончательный ответ

Вместо того чтобы держать 2 массива fields и errorsвы можете держать это в 1 массив и указатель ошибок и значений.

def calc2():
twopi = 2.*np.pi

# annulus of uniform, unit charge density
rmin, rmax = 0.8, 1.2
thmin, thmax = 0, twopi

x = np.arange(-2, 2.1, 0.1)
z = np.arange(-1.55, 1.6, 0.1) # avoid z=0 plane for integration.
X, Z = np.meshgrid(x, z)
s = X.shape

eps = 1E-4
y0 = 0
def Excalc(r, th):
x, y, z = r*np.cos(th), r*np.sin(th), 0.0
return (x0-x) * ((x0-x)**2 + (y0-y)**2 + (z0-z)**2)**-1.5

def Eycalc(r, th):
x, y, z = r*np.cos(th), r*np.sin(th), 0.0
return (y0-y) * ((x0-x)**2 + (y0-y)**2 + (z0-z)**2)**-1.5

def Ezcalc(r, th):
x, y, z = r*np.cos(th), r*np.sin(th), 0.0
return (z0-z) * ((x0-x)**2 + (y0-y)**2 + (z0-z)**2)**-1.5
values = []
for x0, z0 in zip(X.flatten(), Z.flatten()):
# use of lambda in lieu of function calls https://stackoverflow.com/a/49441680/3904031
Ex, Exerr = dblquad(Excalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
Ey, Eyerr = dblquad(Eycalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
Ez, Ezerr = dblquad(Ezcalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
values.append([Ex, Ey, Ez, Exerr, Eyerr, Ezerr])
# errors.append([Exerr, Eyerr, Ezerr])
values = np.array(values).T.reshape(6, *X.shape)

e_vec = values[:3]
e_mag = np.linalg.norm(e_vec, axis=0)

e_vec_err = values[3:]
e_mag_err = np.linalg.norm(e_vec_err, axis=0)

return Evec, Emag, Evecerr, Emagerr

Это привело к очень незначительному снижению времени, так что не избавит вас кучи

lru_cache

рефакторинг calc методы для повторного использования большей части результатом действительно замедлился вещи вниз, это, как я решить, что

from functools import partial, lru_cache

@lru_cache(None)
def calc(r, th, coords_0):
coords_0 = np.array(coords_0)
diff = coords_0 - np.array((r*np.cos(th), r*np.sin(th), 0.0))
return diff * ((diff)**2).sum()**-1.5

def calc3():
rmin, rmax = 0.8, 1.2
thmin, thmax = 0, 2.*np.pi

# s = X.shape

eps = 1E-4
y0 = 0
x = np.arange(-2, 2.1, 0.2)
z = np.arange(-1.55, 1.6, 0.2) # avoid z=0 plane for integration.
X, Z = np.meshgrid(x, z)

calc_axis = lambda r, th, coords_0, axis: calc(r, th, coords_0)[axis]

values = []
for x0, z0 in zip(X.flatten(), Z.flatten()):
coords_0 = (x0, y0, z0)
Excalc, Eycalc, Ezcalc = (partial(calc_axis, coords_0=coords_0, axis=i) for i in range(3))
# use of lambda in lieu of function calls https://stackoverflow.com/a/49441680/3904031
Ex, Exerr = dblquad(Excalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
Ey, Eyerr = dblquad(Eycalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
Ez, Ezerr = dblquad(Ezcalc, thmin, thmax, lambda t: rmin, lambda t: rmax, epsrel=eps)
values.append((Ex, Ey, Ez, Exerr, Eyerr, Ezerr))
print(calc.cache_info())
values = np.array(values).T.reshape(6, *X.shape)

e_vec = values[:3]
e_mag = np.linalg.norm(e_vec, axis=0)

e_vec_err = values[3:]
e_mag_err = np.linalg.norm(e_vec_err, axis=0)
return e_vec, e_mag, e_vec_err, e_mag_err, values

Судя по всему, около половины звонков calc были кэшированы, но движется в расчет этого метода привело к чистому замедление.

numba

Что же в результате увеличения скорости был numba

просто добавив @jit перед def calc(r, th, coords_0): сделан расчет 3 раза быстрее

5
ответ дан 23 марта 2018 в 03:03 Источник Поделиться

Я не уверен, если это так, но если кто-либо из 'калькулятор' функции не называл не раз с тем же r и th значения, тогда functools.lru_cache могли бы помочь ускорить здесь.

Это используется, чтобы украсить метод, и сохраняет его в кэш результаты для каждого набора параметров - возврат кэшированный результат сразу, а не повторного запуска расчета. Он требует больше памяти для хранения результатов, однако, так что вы должны рассмотреть, если что компромисс является приемлемым.

Вы можете быть в состоянии получить некоторую скорость путем слияния основная часть этих методов, а затем вернуть x, y, z и расчетным результатом. Это более вероятно для кэш-памяти, т. к. сейчас у вас 1 кеш, а не 3. Упрощенный пример здесь:

@functools.lru_cache()
def Ecalc(r, th):
x, y, z = r*np.cos(th), r*np.sin(th), 0.0
return x, y, z, ((x0-x)**2 + (y0-y)**2 + (z0-z)**2)**-1.5

# Maybe also lru_cache decorator here?
def new_method(axis, r, th):
x, y, z, result = Ecalc(r, th)
axes = {'x': x, 'y': y, 'z': z}
index = axes[axis]
return (index0 - index) * result

new_method('x', r, th)

1
ответ дан 23 марта 2018 в 09:03 Источник Поделиться