Tutorial: Particle tracking

This example is used to show how to run the particle tracking function. In addition we will compare the times between the 3 different options, 1) Serial in Matlab, 2) Vectorized in Matlab and 3) serial by calling the c++ version of the particle tracking.

For this example we have already prepared the required data for the flow field. The flow domain is a rectangular aquifer with dimensions 40 km x 40 km which has been discretized into 554000 elements and 20 layers

load Data4Part_exmpl1

The data are already prepared in the format suitable for particle tracking. However we need to modify hydraulic conductivity. Note that Kx is a matrix (Nnd x 1), where Nnd is the number of mesh nodes This means that Kx is defined on the nodes. We want to use an anisotropy ratio Kx/Kz = 10 therefore we create the following variable

K{1,1}=Kx;
K{2,1}=Kx/10;
por = 0.25;

where por is the porosity.

For this example we want to generate particle positions so that the particle trajectories are somewhat similar. Therefore we will generate particles from the subdomain which is enclosed by the following coordinates:

partdomX=[2500 7500];
partdomY=[2500 37500];
partdomZ=[-100 -150];

Next we will define the number of particles we want to generate. Note that if you run all the cases this simulation will take few hours because of the large run time of the serial Matlab implementation. You can change the Ncase parameter to specify how many scenarios you want to run (e.g. 4-5 is ok. Choose larger number if you dont need your pc for few hours but not greater than 10)

Nprt=[1 5 10 50 100 250  500 1000 1500 2000];
Ncases = 10;

For the particle tracking we need an option structure, which can be initialized as (while the cursor is on the part_options press F1 to get help about the available options)

part_opt = part_options;

We will keep the default values in this example and change only the mode to 'serial' 'vect' and 'cpp' to run the serial the vectorized and the c++ serial particle tracking. Also we set freqplot a large value to surpress the output.

part_opt.freqplot = 1000;

for k=1:Ncases
    xp=partdomX(1)+(partdomX(2)-partdomX(1))*rand(Nprt(k),1);
    yp=partdomY(1)+(partdomY(2)-partdomY(1))*rand(Nprt(k),1);
    zp=partdomZ(1)+(partdomZ(2)-partdomZ(1))*rand(Nprt(k),1);
    part_opt.mode = 'serial';
    tic;
    [XYZ Vxyz exitflag]=ParticleTracking_main([xp yp zp],XY,Z,TR,TRB,H,K,por,part_opt);
    times(k,1)=toc;
    part_opt.mode = 'vect';
    tic;
    [XYZ Vxyz exitflag]=ParticleTracking_main([xp yp zp],XY,Z,TR,TRB,H,K,por,part_opt);
    times(k,2)=toc;
    part_opt.mode = 'cpp';
    tic;
    [XYZ Vxyz exitflag]=ParticleTracking_main([xp yp zp],XY,Z,TR,TRB,H,K,por,part_opt);
    times(k,3)=toc;
end

Next we plot the time required by each mode to run particle tracking. Obviously the c++ implementation is the fastest. Indeed this is the the option we typically choose. However the other modes can be usefull if one want to get an insight of the particle tracking flow field

plot(Nprt(1:Ncases),times, 'linewidth',1.5)
xlabel('Number of Particles','fontsize',12)
ylabel('Time of particle tracking [sec]','fontsize',12)
hh=legend('Serial','Vectorized','C++');
set(hh,'Location','NorthWest','fontsize',12)

Published with MATLAB® R2012b