sync parallel build functions

Fixes #487

Author: ChrisRackauckas

src/build_function.jl

@@ -140,7 +140,7 @@ function _build_and_inject_function(mod::Module, ex)
     params = typeof(runtimefn).parameters
     fn_expr = GeneralizedGenerated.NGG.from_type(params[3])
 
-    # Inject our externally registered module functions 
+    # Inject our externally registered module functions
     new_expr = ModelingToolkit.inject_registered_module_functions(fn_expr)
 
     # Reconstruct the RuntimeFn's Body
@@ -273,6 +273,7 @@ function _build_function(target::JuliaTarget, rhss, args...;
               end
            end)
         ip_let_expr.args[2] =  ModelingToolkit.build_expr(:block, threaded_exprs)
+		ip_let_expr = :(@sync begin $ip_let_expr end)
     elseif parallel isa DistributedForm
 		numworks = Distributed.nworkers()
 		lens = Int(ceil(length(ip_let_expr.args[2].args)/numworks))
@@ -291,9 +292,11 @@ function _build_function(target::JuliaTarget, rhss, args...;
 		resunpack_exprs = [:($(Symbol(reducevars[iter])) = fetch($(spawnvars[iter]))) for iter in 1:numworks]
 
 		ip_let_expr.args[2] = quote
-			$spawn_exprs
-			$(resunpack_exprs...)
-			$(ip_let_expr.args[2])
+			@sync begin
+				$spawn_exprs
+				$(resunpack_exprs...)
+				$(ip_let_expr.args[2])
+			end
 		end
     elseif parallel isa DaggerForm
         @assert HAS_DAGGER[] "Dagger.jl is not loaded; please do `using Dagger`"
@@ -308,9 +311,11 @@ function _build_function(target::JuliaTarget, rhss, args...;
             $(Symbol(reducevar)) = collect(Dagger.delayed(vcat)($(computevars...)))
         end
         ip_let_expr.args[2] = quote
-            $delayed_exprs
-            $reduce_expr
-            $(ip_let_expr.args[2])
+			@sync begin
+	            $delayed_exprs
+	            $reduce_expr
+	            $(ip_let_expr.args[2])
+			end
         end
 	end
 

test/bigsystem.jl

@@ -46,9 +46,35 @@ end
 f(du,u,nothing,0.0)
 
 multithreadedf = eval(ModelingToolkit.build_function(du,u,parallel=ModelingToolkit.MultithreadedForm())[2])
-_du = rand(N,N,3)
-_u = rand(N,N,3)
-multithreadedf(_du,_u)
+
+# Loop to catch syncronization issues
+for i in 1:100
+   _du = rand(N,N,3)
+   _u = rand(N,N,3)
+   multithreadedf(_du,_u)
+   _du2 = copy(_du)
+   f(_du2,_u,nothing,0.0)
+   @test _du ≈ _du2
+end
+
+#=
+jac = sparse(ModelingToolkit.jacobian(vec(du),vec(u)))
+fjac = eval(ModelingToolkit.build_function(jac,u,parallel=ModelingToolkit.SerialForm())[2])
+multithreadedfjac = eval(ModelingToolkit.build_function(jac,u,parallel=ModelingToolkit.MultithreadedForm())[2])
+
+u = rand(N,N,3)
+J = similar(jac,Float64)
+fjac(J,u)
+
+J2 = similar(jac,Float64)
+multithreadedfjac(J2,u)
+@test J ≈ J2
+
+using FiniteDiff
+J3 = Array(similar(jac,Float64))
+FiniteDiff.finite_difference_jacobian!(J2,(du,u)->f!(du,u,nothing,nothing),u)
+maximum(J2 .- Array(J)) < 1e-5
+=#
 
 using Distributed
 addprocs(4)