Opções para resolver sistemas ODE em GPUs?

14

Eu gostaria de descobrir como resolver sistemas de ODEs em GPUs, em uma configuração 'trivialmente paralelizável'. Por exemplo, fazendo uma análise de sensibilidade com 512 conjuntos de parâmetros diferentes.

Idealmente, eu quero resolver o ODE com um solucionador de tempo inteligente adaptável como CVODE, em vez de um tempo fixo fixo como o Forward Euler, mas executá-lo em uma GPU NVIDIA em vez de CPU.

Alguém já fez isso? Existem bibliotecas para isso?

mirams
fonte
Daí os suportes! Estou pensando em uma técnica baseada na divisão do operador (simulações de eletrofisiologia cardíaca), na qual você resolve EDOs nos nós para obter um termo de origem para o PDE e altera um parâmetro ODE para a próxima iteração.
Mirams
11
É importante especificar se você deseja usar o mesmo intervalo de tempo para cada ODE ou não.
Christian Clason
Além disso, se você estiver especificamente interessado nas equações de bidomain (ou monodomain), poderá dar uma olhada em como o CARP faz isso .
Christian Clason
Timesteps diferentes, se o método for adaptável, ele será necessário para diferentes conjuntos de parâmetros ... obrigado pelo link para o que o CARP está fazendo - um solucionador fixo de timestep Rush Larsen ODE, se eu o ler corretamente.
precisa saber é

Respostas:

6

Você pode procurar na biblioteca odeint do Boost e no Thrust . Eles podem ser combinados conforme discutido aqui .

Juan M. Bello-Rivas
fonte
Isso parece ser um pouco diferente - resolvendo sistemas massivos de ODE na GPU em paralelo (com comunicação). Esse link diz "Experimentamos que o tamanho do vetor paralelo deve ser da ordem de 10 ^ 6 para fazer pleno uso da GPU". Eu estou procurando uma boa maneira de cultivar a O (10) ou O (100) vector dimensionado resolve ODE trivialmente parallelisable ...
Mirams
Você já pensou em escrever diretamente em cuda ou openCL? Se eu entendi direito, o que você está fazendo é iterando sobre alguma equação ODE em cada thread separadamente, não deve ser difícil escrevê-lo diretamente.
Hidro Guy
Eu imagino que seria possível codificar um Forward Euler ou outro método fixo de timestep, em que todo processo de GPU usa o mesmo timestep, com bastante facilidade. Gostaria de saber se alguém conseguiu obter um timestap adaptável como o CVODE funcionando ou se isso é impossível obter eficiência em uma GPGPU?
Mirams
o problema com a gpu é que você precisa escrever um código paralelo a dados. Se você escreve a mesma rotina adaptativa, mas absorve toda essa flexibilidade nos valores de alguns parâmetros, provavelmente é possível codificá-la com eficiência na gpu. Isso também significa que você não pode usar as instruções de ramificação, o que provavelmente é o que você acha que tornaria impossível fazê-lo.
Hydro Guy
11
@mirams, há um exemplo para o odeint que cobre exatamente o que você está procurando: boost.org/doc/libs/1_59_0/libs/numeric/odeint/doc/html/… , veja também github.com/boostorg/odeint/blob/ mestre / exemplos / impulso /… . Além disso, o odeint suporta métodos adaptáveis ​​de várias etapas, como no CVODE: github.com/boostorg/odeint/blob/master/examples/…
Christian Clason
12

A biblioteca DifferentialEquations.jl é uma biblioteca para uma linguagem de alto nível (Julia) que possui ferramentas para transformar automaticamente o sistema ODE em uma versão otimizada para solução paralela em GPUs. Existem duas formas de paralelismo que podem ser empregadas: paralelismo baseado em matriz para sistemas ODE grandes e paralelismo de parâmetros para estudos de parâmetros em sistemas ODE relativamente pequenos (<100). Ele suporta métodos implícitos e explícitos de alta ordem e supera rotineiramente ou combina com outros sistemas em benchmarks (no mínimo, envolve os outros, para que seja fácil verificar e usá-los!)

Para essa funcionalidade específica, você pode dar uma olhada no DiffEqGPU.jl, que é o módulo para o paralelismo automatizado de parâmetros. A biblioteca DifferentialEquations.jl possui funcionalidade para estudos de parâmetros paralelos , e este módulo aumenta as configurações existentes para que o estudo aconteça automaticamente em paralelo. O que se faz é transformar o existente ODEProblem(ou outro DEProblemsemelhante SDEProblem) em um EnsembleProbleme especificar com prob_funccomo os outros problemas são gerados a partir do protótipo. A seguir, são resolvidas 10.000 trajetórias da equação de Lorenz na GPU com um método adaptativo explícito de alta ordem:

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)

Observe que o usuário não precisa escrever nenhum código da GPU e, com um único RTX 2080, esse benchmark é uma melhoria de 5x em relação ao uso de uma máquina Xeon de 16 núcleos com paralelismo multithread. Pode-se, então, verificar o README para saber como executar várias GPUs e executar multiprocessamento + GPUs para utilizar um cluster completo de GPUs simultaneamente . Observe que alternar para multithreading em vez de GPUs é uma alteração de linha: em EnsembleThreads()vez de EnsembleGPUArray().

Então, para solucionadores implícitos, a mesma interface é válida. Por exemplo, o seguinte usa os métodos Rosenbrock de alta ordem e Runge-Kutta implícito:

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)

Embora este formulário exija que você forneça um jacobiano para ser usado na GPU (atualmente, será corrigido em breve), a documentação DifferentialEquations.jl demonstra como fazer cálculos jacobianos simbólicos automáticos em funções numericamente definidas , portanto ainda não há manual trabalhe aqui. Eu recomendo esses algoritmos porque a lógica de ramificação de um método como o CVODE geralmente causa dessincronização de threads e parece não ter um desempenho tão bom quanto o método de Rosenbrock nesses tipos de cenários.

Ao usar DifferentialEquations.jl, você também obtém acesso à biblioteca completa, que inclui funcionalidades como a análise de sensibilidade global, que pode fazer uso dessa aceleração de GPU. Também é compatível com números duplos para análise de sensibilidade local rápida . O código baseado em GPU obtém todos os recursos do DifferentialEquations.jl, como a manipulação de eventos e o grande conjunto de solucionadores de ODE otimizados para diferentes tipos de problemas , o que significa que não é apenas um simples solucionador único de GPU ODE, mas sim um parte de um sistema completo que também possui suporte eficiente à GPU.

Chris Rackauckas
fonte