# Matlab and Julia -- Part II

In the previous post, I explained how to program the solution procedure of Buckley-Leverett equation in Matlab. Here, I'm trying to move everything to Julia. First, you need to install Julia, and a few important packages. Personally, I prefer the last development version. In Ubuntu-based distributions, you can install it by writing the following lines in the terminal.

sudo add-apt-repository ppa:staticfloat/julianightlies sudo add-apt-repository ppa:staticfloat/julia-deps sudo apt-get update sudo apt-get install julia

You can get installers for other OS's here.

First, install a few dependencies, including python 2.7 and matplotlib.
Then, run Julia and add the PyPlot and Roots packages, by typing

Pkg.add("Roots") Pkg.add("PyPlot")

Now, you can proceed to the next step, writing the code in Julia. There are a few things that you need to now when you start working with Julia as a Matlab user. In Matlab, you write

a=linspace(0,1,100); b=1; c=rand(size(a)); a(1) % 0 a(end) % 1 length(a) % 100 size(a) % 1 100 size(b) % 1 1 b*a a.*c

and in Julia, you write

a=linspace(0,1,100) b=1 c=rand(size(a)) a[1] a[end] length(a) # 100 size(a) # (100,) size(b) # () b*a a.*c

In Matlab, all the variables are 2D matrices. In Julia, they are not. Now try to get used to it. For me, it was the hardest part, next to the fact that there is not a good IDE for Julia yet.

No I convert everything right away to Julia, and you can compare it to the Matlab code for yourself.

# A simple compact Julia script to solve Buckley-Leverett equation # for water flooding. # Written by AAE, TU Delft, September 2014 # Any sort of reproduction and redistribution of this code is highly encouraged. # Note: The extrapolation of relperm curves is not included to keep # the code as simple as possible. # run the following two lines first: # using Roots # using PyPlot # viscosities muw = 1e-3 muo = 10e-3 # initialize the Corey model parameters kro0 = 1 krw0 = 1 no = 2 nw = 4 sor = 0.1 swc = 0.2 u = 1e-3 phi = 0.3 # define the functions sws(sw)=((sw-swc)/(1-swc-sor)) dsws = 1/(1-swc-sor) krw(sw)=(krw0*sws(sw).^nw) kro(sw)=(kro0*(1-sws(sw)).^no) # derivatives dkrw(sw)=(dsws*nw*krw0*sws(sw).^(nw-1)) dkro(sw)=(-dsws*no*kro0*(1-sws(sw)).^(no-1)) fw(sw)=((krw(sw)/muw)./(krw(sw)/muw+kro(sw)/muo)) dfw(sw)=((dkrw(sw)/muw.*kro(sw)/muo-krw(sw)/muw.*dkro(sw)/muo)./(krw(sw)/muw+kro(sw)/muo).^2) # define initial conditions sw0 = swc # the reservoir is filled with oil and connate water # solve the nl equation to find the shock front saturation f_shock(sw)=(dfw(sw)-(fw(sw0)-fw(sw))./(sw0-sw)) sw_shock = fzero(f_shock, 0.5) # injection saturation sw_inj = 1-sor # not to make things too complicated # water saturation between injection and shock sw_rf = linspace(sw_inj, sw_shock, 100) sw_all = linspace(swc, 1-sor, 100) # whole range of possible sw xt_all = u/phi*dfw(sw_all) # solution over a whole possible range of sw xt_rf=u/phi*dfw(sw_rf) xt_shock = xt_rf[end] figure(1) subplot(3,1,1) plot(sw_all, krw(sw_all), sw_all, kro(sw_all)) axis([0, 1, 0, 1]) xlabel(L"S_w") ylabel("Relative permeabilities") legend((L"k_{rw}", L"k_{ro}")) subplot(3,1,2) plot(sw_all, fw(sw_all)) axis([0, 1, 0, 1]) xlabel(L"S_w") ylabel(L"F_w") subplot(3,1,3) plot([xt_rf, xt_shock, maximum(xt_all)], [sw_rf, sw0, sw0], "r", xt_all, sw_all, "--b") xlabel("x/t") ylabel(L"S_w") legend(("BL solution",))

Now is the important step, running the code. Load the packages that you are going to call from your code:

using Roots using PyPlot

We will use `fzero`

function from `Roots`

package to find the shock front. We will also use the `PyPlot`

package to visualize the results. Now load the code in `Julia`

by typing

reload("BLjulia.jl")

If everything is fine, it will show you the following plot:

I'll discuss some of the small issues that I had with Julia as a Matlab user.

## Comments

Comments powered by Disqus