Варианты решения систем ODE на графических процессорах?

15

Я хотел бы вывести системы решения ODE на графические процессоры в «тривиально распараллеливаемой» обстановке. Например, анализ чувствительности с 512 различными наборами параметров.

В идеале я хочу решать ODE с помощью интеллектуального адаптивного временного шага, такого как CVODE, а не с фиксированным временным шагом, как Forward Euler, но запускать его на GPU NVIDIA вместо CPU.

Кто-нибудь делал это? Есть библиотеки для этого?

mirams
источник
Отсюда и скобки! Я рассматриваю метод, основанный на разделении операторов (моделирование электрофизиологии сердца), в котором вы решаете ODE в узлах, чтобы получить исходный термин для PDE, а затем меняете параметр ODE для следующей итерации.
Мирамс
1
Важно указать, хотите ли вы использовать один и тот же временной шаг для каждого ODE или нет.
Кристиан Клэйсон
Кроме того, если вы особенно заинтересованы в уравнениях бидомена (или монодомена), вы можете посмотреть, как это делает CARP .
Кристиан Клэйсон,
Разные временные шаги, если метод адаптивный, тогда он понадобится для разных наборов параметров ... спасибо за ссылку на то, что делает CARP - фиксированный временный шаг для решения ODE Rush Larsen, если я правильно его прочитал.
Мирамс

Ответы:

6

Вы можете заглянуть в odeint библиотеку Boost и Thrust . Они могут быть объединены как обсуждено здесь .

Хуан М. Белло-Ривас
источник
Кажется, что это немного по-другому - решение больших систем ODE на GPU параллельно (с обменом данными). Эта ссылка гласит: «Мы убедились, что размер вектора, над которым распараллеливаются, должен быть порядка 10 ^ 6, чтобы в полной мере использовать графический процессор». Я ищу хороший способ выработки O (10) или O (100) векторно-размерных тривиально распараллеливающихся
решений
Вы думали о написании непосредственно в CUDA или OpenCL? Если я правильно понял, то, что вы делаете - это перебираете какое-то уравнение ODE в каждом потоке в отдельности, не должно быть сложным написать его напрямую.
Hydro Guy
Я полагаю, что можно было бы кодировать Forward Euler или другой метод с фиксированным временным шагом, где каждый процесс GPU довольно легко использует один и тот же временной шаг, я хотел бы знать, удалось ли кому-нибудь получить адаптивный временной шаг, такой как CVODE, или это невозможно сделать работоспособным на GPGPU?
Мирамс
проблема с gpu в том, что вам нужно писать параллельный к данным код. Если вы пишете ту же адаптивную подпрограмму, но впитываете всю эту гибкость в значения некоторых параметров, возможно, ее можно эффективно закодировать в gpu. Это также означает, что вы не можете использовать ветвление в инструкциях, что, вероятно, и делает это невозможным.
Hydro Guy
1
@mirams есть пример для odeint, который охватывает именно то, что вы ищете: boost.org/doc/libs/1_59_0/libs/numeric/odeint/doc/html/… , см. также github.com/boostorg/odeint/blob/ мастер / примеры / тяга /… . Кроме того, odeint поддерживает адаптивные многошаговые методы, как в CVODE: github.com/boostorg/odeint/blob/master/examples/…
Кристиан Клэйсон,
12

Библиотека DifferentialEquations.jl - это библиотека для языка высокого уровня (Julia), которая имеет инструменты для автоматического преобразования системы ODE в оптимизированную версию для параллельного решения на графических процессорах. Можно использовать две формы параллелизма: параллелизм на основе массива для больших систем ODE и параллелизм параметров для исследования параметров на относительно небольших (<100) системах ODE. Он поддерживает неявные и явные методы высокого порядка и обычно превосходит или сопоставляет другие системы в тестах (по крайней мере, он оборачивает другие, поэтому их легко проверить и использовать!)

Для этой специфической функциональности вы можете взглянуть на DiffEqGPU.jl, который является модулем для автоматического параллелизма параметров. Библиотека DifferentialEquations.jl имеет функциональность для параллельных исследований параметров , и этот модуль дополняет существующие конфигурации для автоматического параллельного изучения. То, что вы делаете, это преобразовывает их существующие ODEProblem(или другие DEProblemподобные SDEProblem) в EnsembleProblemи указывает, prob_funcкак другие проблемы генерируются из прототипа. Следующее решает 10 000 траекторий уравнения Лоренца на графическом процессоре с помощью явного адаптивного метода высокого порядка:

using OrdinaryDiffEq, DiffEqGPU
function lorenz(du,u,p,t)
 @inbounds begin
     du[1] = p[1]*(u[2]-u[1])
     du[2] = u[1]*(p[2]-u[3]) - u[2]
     du[3] = u[1]*u[2] - p[3]*u[3]
 end
 nothing
end

u0 = Float32[1.0;0.0;0.0]
tspan = (0.0f0,100.0f0)
p = (10.0f0,28.0f0,8/3f0)
prob = ODEProblem(lorenz,u0,tspan,p)
prob_func = (prob,i,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func)
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10_000,saveat=1.0f0)

Обратите внимание, что пользователю не нужно писать код графического процессора, и с одним RTX 2080 этот тест оценивается как пятикратное улучшение по сравнению с использованием 16-ядерного компьютера Xeon с многопоточным параллелизмом. Затем можно проверить README о том, как использовать несколько графических процессоров и мультипроцессорные + графические процессоры для одновременного использования полного кластера графических процессоров . Обратите внимание, что переключение на многопоточность вместо графических процессоров является одним изменением строки: EnsembleThreads()вместоEnsembleGPUArray() .

Тогда для неявных решателей имеет место тот же интерфейс. Например, следующее использует Rosenbrock высокого порядка и неявные методы Runge-Kutta:

function lorenz_jac(J,u,p,t)
 @inbounds begin
     σ = p[1]
     ρ = p[2]
     β = p[3]
     x = u[1]
     y = u[2]
     z = u[3]
     J[1,1] = -σ
     J[2,1] = ρ - z
     J[3,1] = y
     J[1,2] = σ
     J[2,2] = -1
     J[3,2] = x
     J[1,3] = 0
     J[2,3] = -x
     J[3,3] = -β
 end
 nothing
end

function lorenz_tgrad(J,u,p,t)
 nothing
end

func = ODEFunction(lorenz,jac=lorenz_jac,tgrad=lorenz_tgrad)
prob_jac = ODEProblem(func,u0,tspan,p)
monteprob_jac = EnsembleProblem(prob_jac, prob_func = prob_func)

@time solve(monteprob_jac,Rodas5(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10_000,saveat=1.0f0)
@time solve(monteprob_jac,TRBDF2(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10_000,saveat=1.0f0)

Хотя эта форма требует, чтобы вы указали якобиан для использования в графическом процессоре (в настоящее время это будет исправлено в ближайшее время), в документации DifrentialEquations.jl показано, как выполнять автоматические символьные вычисления якобиана для численно определенных функций. , поэтому руководство по-прежнему отсутствует. работа здесь. Я бы настоятельно рекомендовал эти алгоритмы, потому что логика ветвления такого метода, как CVODE, обычно вызывает рассинхронизацию потоков и, похоже, в любом случае не работает так же хорошо, как метод Розенброка в этих типах сценариев.

Используя diffrentialEquations.jl, вы также получаете доступ ко всей библиотеке, которая включает в себя такие функции, как глобальный анализ чувствительности, который может использовать это ускорение GPU. Он также совместим с двойными числами для быстрого анализа локальной чувствительности . Код на основе графического процессора получает все функции DifrentialEquations.jl, такие как обработка событий и большой набор решателей ODE, которые оптимизированы для различных типов проблем , то есть это не просто одноразовый GPU ODE решатель, а вместо этого часть полнофункциональной системы, которая также имеет эффективную поддержку графического процессора.

Крис Ракауцкас
источник