Tutorial on Parallelization

This tutorial shows how to use jInv to parallelize evaluations of your forward operator. We will see three modes for coputing the forward problems

Putting all sources into a single ForwardProbType gives a single forward problem that is evaluated by one worker.
Creating multiple instances of ForwardProbType each containing a subset of sources and storing them in an Array the forward problems are solved in parallel and the assignment of problems to workers is created dynamically.
Creating one ForwardProbType per worker and dividing the sources among them yields an array of RemoteChannel objects. The assignment is created prior to running the problem and then kept fixed. We refer to this as static scheduling.

All of the above strategies give the same results and can be used with the same syntax getData(model,pFor) thanks to Julia's multiple dispatch. The optimal choice depends on the characteristics of the problem at hand as well as on the available computational resources.

References

Please refer to [1] for a detailed description of the parallelization strategies in jInv. The fast marching method used to solve the Eikonal equation is presented in [2]. The test data used here is generated using jInv and based on the 3D SEG/EAGE model of a salt reservoir described in [3].

[1] Ruthotto L, Treister E, Haber E: jInv - a flexible Julia package for PDE parameter estimation, arXiv:1606.0739 [cs.MS] 2016

[2] Treister E, Haber E: A fast marching algorithm for the factored eikonal equation, Journal of Computational Physics, 324(1), p. 210–225, 2016

[3] F. Aminzadeh, B. Jean, and T. Kunz. 3-D salt and overthrust models. Society of Exploration Geophysicists, 1997.

Add workers and load required packages

Before we start, we should add workers to the session. Here, we add local workers, however, the workers added could also be on remote computers (see also Julia Documentation)



In [1]:

    
addprocs(8)









    Out[1]:





8-element Array{Int64,1}:
 2
 3
 4
 5
 6
 7
 8
 9



In [2]:

    
using MAT
using EikonalInv
using jInv.Mesh
using jInv.ForwardShare
using jInv.Utils
using Test

Load test data and setup mesh

Here, we use the test data for a travel time tomography experiment stored in exEikonal.mat



In [3]:

    
matfile   = matread("exEikonal.mat")
Receivers = matfile["Receivers"]# sparse matrix of discrete receivers
Sources   = matfile["Sources"] # sparse matrix of discrete sources
slowref   = matfile["slowref"] # reference slowness, serves as starting guess;
domain    = matfile["domain"]
n         = matfile["n"]

# set up mesh
h           = (domain[2:2:end]-domain[1:2:end])./n
domainNodal = copy(domain); domainNodal[1:2:end] += h./2; domainNodal[2:2:end] -= h./2
Mfwd        = getRegularMesh(domain,n-1);

Option 1) Single Forward Problem -> Serial Computation

The simplest way to evaluate the forward operator and simulate data is by generating a single EikonalInvParam containing all sources. This way, there is no parallelization except the one included in the method for solving the forward problem, which in the case of fast marching means almost no parallelization.



In [5]:

    
pFor = getEikonalInvParam(Mfwd,Sources,Receivers,false,false);

@time dobs,pFor = getData(slowref,pFor);









    



  7.562552 seconds (97.03 M allocations: 1.527 GB, 2.13% gc time)

Option 2) Array of Forward Problems --> Dynamic Scheduling

The second, more advanced option is to divide the sources into different instances of EikonalInvParam and use dynamic scheduling. Note that simulating the data can be done independently for different sources. Thus, when computing the data, the problems are spawned out dynamically to the next available worker. Thereby latency times are limited, however, all data required to solve the problem (including meshes, differential operators, etc) need to be sent to the remote worker.



In [7]:

    
pForp = Array{EikonalInvParam}(size(Sources,2))

for k=1:size(Sources,2)
    pForp[k] = getEikonalInvParam(Mfwd,SparseMatrixCSC(Sources[:,k]),Receivers,false,false)
end
@time begin
    # solve forward problems in parallel
    dobspFuture,pForp = getData(slowref,pForp);

    # fetch results from workers
    dobsp = zeros(size(Receivers,2),size(Sources,2))
    for k=1:size(Sources,2)
        dobsp[:,k] = fetch(dobspFuture[k])
    end
end

# make sure that results corresponds to serial computation
@test norm(dobsp-dobs)/norm(dobs) < 1e-14









    



  2.411892 seconds (56.13 k allocations: 86.586 MB, 1.15% gc time)






    Out[7]:





Test Passed
  Expression: norm(dobsp - dobs) / norm(dobs) < 1.0e-14
   Evaluated: 0.0 < 1.0e-14

Option 3) Array of RemoteChannels --> Static Scheduling

A third option is to divide the sources among the available workers before calling the forward problem and then keeping the assignmnent fixed. We refer to this as static scheduling. Here, the EikonalInvParam is actually initialized on the remote worker and does not need to be sent back and forth. This way communication is kept at a minimum, but latency times might occur if some problems take longer to solve than others.



In [9]:

    
pForr = Array{RemoteChannel}(nworkers())


# create forward problems on remote workers
SourcesSubInd = Array{Any}(nworkers())
for k=1:nworkers()
    SourcesSubInd[k] = collect(k:nworkers():size(Sources,2))
    pForr[k] = initRemoteChannel(getEikonalInvParam, workers()[k], Mfwd,Sources[:,SourcesSubInd[k]],Receivers,false,false) 
end

@time begin
    # solve forward problems in parallel
    dobsrFuture,pForr = getData(slowref,pForr);
    
    # fetch results
    dobsr = zeros(size(Receivers,2),size(Sources,2))
    for k=1:nworkers()
        dobsr[:,SourcesSubInd[k]] = fetch(dobsrFuture[k])
    end   
end

# make sure that results corresponds to serial computation
@test norm(dobsr-dobs)/norm(dobs) < 1e-14









    



  2.373038 seconds (3.26 k allocations: 2.151 MB)






    Out[9]:





Test Passed
  Expression: norm(dobsr - dobs) / norm(dobs) < 1.0e-14
   Evaluated: 0.0 < 1.0e-14



In [ ]: