Merge pull request #40 from michakraus/dev-parallel

Fix some parallel simulation issues

Project.toml

@@ -1,7 +1,7 @@
 name = "GeometricIntegrators"
 uuid = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
 authors = ["Michael Kraus <[email protected]>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 DecFP = "55939f99-70c6-5e9b-8bb0-5071ed7d61fd"

prototyping/profile_simulation.jl

@@ -0,0 +1,49 @@
+# run with julia --track-allocation=user
+
+using GeometricIntegrators
+using GeometricIntegrators.TestProblems.HarmonicOscillatorProblem
+using Profile
+using Test
+
+set_config(:nls_atol_break, Inf)
+set_config(:nls_rtol_break, Inf)
+set_config(:nls_stol_break, Inf)
+
+
+Δt = 0.1
+nt = 10000
+ns = 100
+
+nsave = 1
+nwrte = 10000
+
+h5file = "test.hdf5"
+
+ode = harmonic_oscillator_ode(vcat(rand(1,ns), zeros(1,ns)))
+#tab = getTableauImplicitMidpoint()
+tab = getTableauERK4()
+
+
+sim = Simulation(ode, tab, Δt, "Serial Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)
+run!(sim)
+rm(h5file)
+
+sim = Simulation(ode, tab, Δt, "Serial Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)
+
+Profile.clear()
+Profile.clear_malloc_data()
+
+@profile run!(sim)
+rm(h5file)
+
+sim = ParallelSimulation(ode, tab, Δt, "Parallel Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)
+run!(sim)
+rm(h5file)
+
+sim = ParallelSimulation(ode, tab, Δt, "Parallel Harmonic Oscillator Test", h5file, nt; nsave=nsave, nwrite=nwrte)
+
+Profile.clear()
+Profile.clear_malloc_data()
+
+@profile run!(sim)
+rm(h5file)

src/CommonFunctions.jl

@@ -22,10 +22,14 @@ module CommonFunctions
 
     write_to_hdf5() = nothing
 
-
     export reset!, cut_periodic_solution!
 
     reset!() = nothing
     cut_periodic_solution!() = nothing
 
+    export eachsample, eachtimestep
+
+    eachsample() = nothing
+    eachtimestep() = nothing
+
 end

src/Solutions.jl

@@ -39,7 +39,6 @@ module Solutions
 
     export Solution, ParallelSolution, DeterministicSolution, StochasticSolution
     export nsave, nsamples, ntime, timesteps, offset, conv, hdf5
-    export eachsample, eachtimestep
     export get_solution, get_solution!, set_solution!
 
     include("solutions/solution.jl")

src/integrators/rk/integrators_dirk.jl

@@ -36,10 +36,11 @@ mutable struct ParametersDIRK{DT, TT, ET <: ODE{DT,TT}, D, S} <: Parameters{DT,T
 
     t::TT
     q::Vector{DT}
+    V::Vector{Vector{DT}}
 end
 
 function ParametersDIRK(equ::ET, tab::TableauDIRK{TT}, Δt::TT) where {DT, TT, ET <: ODE{DT,TT}}
-    ParametersDIRK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, 0, zeros(DT, equ.d))
+    ParametersDIRK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, 0, zeros(DT, equ.d), create_internal_stage_vector(DT, equ.d, tab.s))
 end
 
 
@@ -143,47 +144,45 @@ function initial_guess!(int::IntegratorDIRK, sol::AtomicSolutionODE)
 end
 
 
-function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
+function compute_stages!(x::Vector{ST}, Q::Vector{ST}, V::Vector{ST}, Y::Vector{ST},
                                     params::ParametersDIRK{DT,TT,ET,D,S}, i::Int) where {ST,DT,TT,ET,D,S}
 
     local tᵢ::TT
 
-    @assert S == length(Q) == length(V) == length(Y)
-    @assert D == length(Q[1]) == length(V[1]) == length(Y[1])
+    @assert D == length(Q) == length(V) == length(Y)
 
     # copy x to Y and compute Q = q + Δt Y
     for k in 1:D
-        Y[i][k] = x[k]
-        Q[i][k] = params.q[k] + params.Δt * Y[i][k]
+        Y[k] = x[k]
+        Q[k] = params.q[k] + params.Δt * Y[k]
     end
 
     # compute V = v(Q)
     tᵢ = params.t + params.Δt * params.tab.q.c[i]
-    params.equ.v(tᵢ, Q[i], V[i])
+    params.equ.v(tᵢ, Q, V)
 end
 
 
 "Compute stages of fully implicit Runge-Kutta methods."
-@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersDIRK{DT,TT,ET,D,S}, i::Int) where {ST,DT,TT,ET,D,S}
-
-    cache = IntegratorCacheDIRK{ST,D,S}()
+function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersDIRK{DT,TT,ET,D,S}, i::Int) where {ST,DT,TT,ET,D,S}
+    local Q = zeros(ST,D)
+    local V = zeros(ST,D)
+    local Y = zeros(ST,D)
 
-    quote
-        compute_stages!(x, $cache.Q, $cache.V, $cache.Y, params, i)
+    local y1::ST
+    local y2::ST
 
-        local y1::ST
-        local y2::ST
+    compute_stages!(x, Q, V, Y, params, i)
 
-        # compute b = - (Y-AV)
-        for k in 1:D
-            y1 = 0
-            y2 = 0
-            for j in 1:S
-                y1 += params.tab.q.a[i,j] * $cache.V[j][k]
-                y2 += params.tab.q.â[i,j] * $cache.V[j][k]
-            end
-            b[k] = - $cache.Y[i][k] + (y1 + y2)
+    # compute b = - (Y-AV)
+    for k in 1:D
+        y1 = params.tab.q.a[i,i] * V[k]
+        y2 = params.tab.q.â[i,i] * V[k]
+        for j in 1:i-1
+            y1 += params.tab.q.a[i,j] * params.V[j][k]
+            y2 += params.tab.q.â[i,j] * params.V[j][k]
         end
+        b[k] = - Y[k] + (y1 + y2)
     end
 end
 
@@ -211,7 +210,10 @@ function integrate_step!(int::IntegratorDIRK{DT,TT}, sol::AtomicSolutionODE{DT,T
         check_solver_status(int.solver[i].status, int.solver[i].params)
 
         # compute vector field at internal stages
-        compute_stages!(int.solver[i].x, int.cache.Q, int.cache.V, int.cache.Y, int.params, i)
+        compute_stages!(int.solver[i].x, int.cache.Q[i], int.cache.V[i], int.cache.Y[i], int.params, i)
+
+        # copy velocity field to parameters
+        int.params.V[i] .= int.cache.V[i]
     end
 
     # compute final update

src/integrators/rk/integrators_firk.jl

@@ -189,27 +189,29 @@ function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{
 end
 
 "Compute stages of fully implicit Runge-Kutta methods."
-@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFIRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
+function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFIRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
+    # temporary variables
+    local y1::ST
+    local y2::ST
 
-    cache = IntegratorCacheFIRK{ST,D,S}()
+    # create caches for internal stages
+    Q = create_internal_stage_vector(ST, D, S)
+    V = create_internal_stage_vector(ST, D, S)
+    Y = create_internal_stage_vector(ST, D, S)
 
-    quote
-        compute_stages!(x, $cache.Q, $cache.V, $cache.Y, params)
+    # compute stages from nonlinear solver solution x
+    compute_stages!(x, Q, V, Y, params)
 
-        local y1::ST
-        local y2::ST
-
-        # compute b = - (Y-AV)
-        for i in 1:S
-            for k in 1:D
-                y1 = 0
-                y2 = 0
-                for j in 1:S
-                    y1 += params.tab.q.a[i,j] * $cache.V[j][k]
-                    y2 += params.tab.q.â[i,j] * $cache.V[j][k]
-                end
-                b[D*(i-1)+k] = - $cache.Y[i][k] + params.Δt * (y1 + y2)
+    # compute b = - (Y-AV)
+    for i in 1:S
+        for k in 1:D
+            y1 = 0
+            y2 = 0
+            for j in 1:S
+                y1 += params.tab.q.a[i,j] * V[j][k]
+                y2 += params.tab.q.â[i,j] * V[j][k]
             end
+            b[D*(i-1)+k] = - Y[i][k] + params.Δt * (y1 + y2)
         end
     end
 end

src/integrators/rk/integrators_iprk.jl

@@ -153,50 +153,56 @@ end
 
 
 "Compute stages of implicit partitioned Runge-Kutta methods."
-@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersIPRK{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}
-    cache = IntegratorCacheIPRK{ST,D,S}()
-
-    quote
-        compute_stages!(x, $cache, params)
-
-        # compute b = - [(Y-AV), (Z-AF)]
-        for i in 1:S
-            for k in 1:D
-                b[2*(D*(i-1)+k-1)+1] = - $cache.Y[i][k]
-                b[2*(D*(i-1)+k-1)+2] = - $cache.Z[i][k]
-                for j in 1:S
-                    b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * $cache.V[j][k]
-                    b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * $cache.F[j][k]
-                end
+function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersIPRK{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}
+    # create caches for internal stages
+    Q = create_internal_stage_vector(ST, D, S)
+    P = create_internal_stage_vector(ST, D, S)
+    V = create_internal_stage_vector(ST, D, S)
+    F = create_internal_stage_vector(ST, D, S)
+    Y = create_internal_stage_vector(ST, D, S)
+    Z = create_internal_stage_vector(ST, D, S)
+
+    # compute stages from nonlinear solver solution x
+    compute_stages!(x, Q, V, Y, P, F, Z, params)
+
+    # compute b = - [(Y-AV), (Z-AF)]
+    for i in 1:S
+        for k in 1:D
+            b[2*(D*(i-1)+k-1)+1] = - Y[i][k]
+            b[2*(D*(i-1)+k-1)+2] = - Z[i][k]
+            for j in 1:S
+                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * V[j][k]
+                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * F[j][k]
             end
         end
     end
 end
 
 
-function compute_stages!(x::Vector{ST}, cache::IntegratorCacheIPRK{ST,D,S},
+function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
+                                        P::Vector{Vector{ST}}, F::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                         params::ParametersIPRK{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}
     local tqᵢ::TT
     local tpᵢ::TT
 
     for i in 1:S
         for k in 1:D
             # copy y to Y and Z
-            cache.Y[i][k] = x[2*(D*(i-1)+k-1)+1]
-            cache.Z[i][k] = x[2*(D*(i-1)+k-1)+2]
+            Y[i][k] = x[2*(D*(i-1)+k-1)+1]
+            Z[i][k] = x[2*(D*(i-1)+k-1)+2]
 
             # compute Q and P
-            cache.Q[i][k] = params.q[k] + params.Δt * cache.Y[i][k]
-            cache.P[i][k] = params.p[k] + params.Δt * cache.Z[i][k]
+            Q[i][k] = params.q[k] + params.Δt * Y[i][k]
+            P[i][k] = params.p[k] + params.Δt * Z[i][k]
         end
 
         # compute time of internal stage
         tqᵢ = params.t + params.Δt * params.tab.q.c[i]
         tpᵢ = params.t + params.Δt * params.tab.p.c[i]
 
         # compute v(Q,P) and f(Q,P)
-        params.equ.v(tqᵢ, cache.Q[i], cache.P[i], cache.V[i])
-        params.equ.f(tpᵢ, cache.Q[i], cache.P[i], cache.F[i])
+        params.equ.v(tqᵢ, Q[i], P[i], V[i])
+        params.equ.f(tpᵢ, Q[i], P[i], F[i])
     end
 end
 
@@ -259,7 +265,8 @@ function integrate_step!(int::IntegratorIPRK{DT,TT}, sol::AtomicSolutionPODE{DT,
     check_solver_status(int.solver.status, int.solver.params)
 
     # compute vector fields at internal stages
-    compute_stages!(int.solver.x, int.cache, int.params)
+    compute_stages!(int.solver.x, int.cache.Q, int.cache.V, int.cache.Y,
+                                  int.cache.P, int.cache.F, int.cache.Z, int.params)
 
     # compute final update
     update_solution!(sol.q, sol.q̃, int.cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))

src/problems/harmonic_oscillator.jl

@@ -6,11 +6,11 @@ module HarmonicOscillatorProblem
            harmonic_oscillator_dae, harmonic_oscillator_idae, harmonic_oscillator_pdae
 
 
-    Δt  = 0.1
-    nt  = 10
+    const Δt = 0.1
+    const nt = 10
 
-    k = 0.5
-    ω = √k
+    const k = 0.5
+    const ω = √k
 
     t₀=0.0
     q₀=[0.5, 0.0]