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
```