Solving ODE with inheritance of initial condition with parallelization by threads

General Usage
1

Hello!
I have an ODE and I need to perform inheritance of initial conditions when changing two parameters. I can do this without parallelizing threads, and everything works correctly. But when I start parallelizing, I get the wrong result. The variable u0_start_2d takes incorrect values. I assume that one thread interferes with the work of another thread. But I don’t understand how this happens.

A fragment of code that I have paralleled:

Threads.@threads for p2_index_range in eachindex(p2_range)
        for p1_index_range in eachindex(p1_range)

            if p1_index_range == 1
                u0_start_2d = SVector{6}(matrix_u0_end[p1_index_range, p2_index_range, :]);
                continue;
            end
            
            loc_params = copy(params);
            loc_params[index_p1] = p1_range[p1_index_range];
            loc_params[index_p2] = p2_range[p2_index_range];

            prob = ODEProblem(sys, u0_start_2d, tspan, loc_params);
            sol = solver(prob, integ_set);

            if sol.retcode == ReturnCode.MaxIters
                exitcode = 1;
            end

            u0_end_2d = sol[end];

            matrix_u0_start[p1_index_range, p2_index_range, :] = u0_start_2d;
            matrix_u0_end[p1_index_range, p2_index_range, :] = u0_end_2d;
            matrix_params[p1_index_range, p2_index_range, :] = [loc_params[index_p1], loc_params[index_p2]]
            
            u0_start_2d = u0_end_2d;

            if exitcode == 1
                exit();
            end;
        end
    end

Link on full project: source code - Google Drive

2

So something from before the loop might be interesting but I guess you do something like:

u0_start_2d = ....
Threads.@threads for p2_index_range in eachindex(p2_range)
        for p1_index_range in eachindex(p1_range)

            if p1_index_range == 1
                u0_start_2d = SVector{6}(matrix_u0_end[p1_index_range, p2_index_range, :]);
                continue;
            end
            ...
            u0_start_2d = u0_end_2d;

Which means each thread sees the same variable u0_start_2d and overwrite each others values!
A fix can simply be moving the variable inside the outer loop:

### u0_start_2d = .... # from here
Threads.@threads for p2_index_range in eachindex(p2_range)
        u0_start_2d = .... # to here
        for p1_index_range in eachindex(p1_range)

then the variable will be local in each thread and is safe to use and overwrite.

You could also refactor this a bit to:

Threads.@threads for p2_index_range in eachindex(p2_range)
        u0_start_2d = SVector{6}(matrix_u0_end[1, p2_index_range, :])
        for p1_index_range in eachindex(p1_range)
           p1_index_range == 1 && continue

           loc_params = copy(params)
           # rest of code...

Also note you don’t need “;” at the end of each line - we are not in C/Java land anymore and thus free of the oppression of semicolons :smile:

1 Like
3

Thank you very much
Everything is working now

4

@Sergey_Novak If you are not aware of if, you may be interested in the EnsembleProblem interface. https://docs.sciml.ai/DiffEqDocs/dev/features/ensemble/