Values of arrays changed unexpectedly due to storage allocation issue?

General Usage
1

Not sure whether I am asking this question in the correct section.
I wanted to build a solver for the 1D advection-dispersion-equation, with Dirichlet boundary condition at the influx end and Neumann BC at the efflux end. The code below was just for a test so there are some hardcoding.
First, a domain with n nodes was initiated as A. Then the function SolveNADE() solves the ADE numerically once, by delta Δt. As I have expected, the initial array A should not be changed since I am traversing the array C_over_Domain. However, this is not the case. The A would have new values after each for i in 1:5000000 loop. What could be the reason and how should I avoid this? Thanks.

nodesN = 100
A = zeros(nodesN)

# the physical condition, alias: PhyCond in functions
struct PhysicalCondition
    u::Float64  # advection velocity [m/s]
    D::Float64  # diffusion coefficient [m²/s]
end
# numerical scheme, alias: Scheme in functions
struct NumericalScheme
    Δx::Float64
    Δt::Float64
end

𝒖 = 4e-6 # advection velocity [m/s]
𝑫 = 1e-9 # diffusion coefficient [m²/s]
phycondition = PhysicalCondition(𝒖, 𝑫)
dx = 1
dt = 1
numericscheme = NumericalScheme(dx, dt)

function SolveNADE(PhyCond, Scheme, domain)
    C_next = domain
    C_next[1] = 1
    for i in (2:size(C_next)[1]-1)
        C_next[i] = (- PhyCond.u * Scheme.Δt / Scheme.Δx * (domain[i] - domain[i-1]) 
                        + PhyCond.D * Scheme.Δt / Scheme.Δx^2 * (domain[i-1] - 2*domain[i] + domain[i+1])
                        + domain[i])
    end
    C_next[end] = (- PhyCond.u * Scheme.Δt / Scheme.Δx * (domain[end] - domain[end-1]) 
                        + PhyCond.D * Scheme.Δt / Scheme.Δx^2 * (domain[end-1] - domain[end])
                        + domain[end])
    return C_next
end

C_over_Domain = A
for i in 1:5000000
    new_domain = C_over_Domain
    C_over_Domain = SolveNADE(phycondition, numericscheme, new_domain)
end
2

You should probably copy() here:

C_over_Domain = copy(A)
2 Likes
3

Thanks! Now I´ve got the point.

4

Hope these further discussions & documentations could help other newcomers of Julia like me.

https://docs.julialang.org/en/v1/base/base/#=

5

Since a (not-so) old question of mine was tagged, this is related: Assignment and mutation · JuliaNotes.jl

1 Like
6

A further question regarding “Assignment and Mutation”:
Why is this happening? To obtain the results that I´m expecting, for each loop I have to make a copy.

x = [1, 3, 6]
emptysaver = []
for i in 1:3
    x[i] = x[i] + 10
    push!(emptysaver, x)
end
emptysaver

3-element Vector{Any}:
 [11, 13, 16]
 [11, 13, 16]
 [11, 13, 16]

x = [1, 3, 6]
emptysaver = []
for i in 1:3
    x[i] = x[i] + 10
    new = copy(x)
    push!(emptysaver, new)
end
emptysaver

3-element Vector{Any}:
 [11, 3, 6]
 [11, 13, 6]
 [11, 13, 16]
7

Vectors are passed by sharing, so when you push a vector to a collection of vectors, you are passing the reference to that vector. Thus, the reference to the same vector is being pushed every time.

To push the reference of a new vector, you need to first create the new vector, by copying.

1 Like
9

It is often helpful to write this as a comprehension, which will create new object (rather than copy references) and also yield the right type (note your vector is a Vector{Any} which kills performance, don’t initialise vectors with []):

julia> [x .+ 10 for _ ∈ 1:3]
3-element Vector{Vector{Int64}}:
 [11, 13, 16]
 [11, 13, 16]
 [11, 13, 16]
2 Likes