# Numerical-Code: Warnings and Errors due to scopes

I am a newbie in Julia an I have written the following code:

``````using LinearAlgebra
using Plots

t = 0;

epsilon = 0.1;

########## initial data: Sod problem (rho,u,T) ##########

Prim_l = [1 0 1];                           # 0 <= x <= 0.5
Prim_r = [0.125 0 0.1];                     # 0.5 < x <= 1

########## domain in space ##########

x_start = 0;                                # start
x_end   = 1;                                # end

k   = 250;                                  # cell discretization, grid-points
dx  = (x_end-x_start)/k;

########## velocity space (v = molecular velocities at certain point x at time t) ##########

v_start = -4.5;                             # start
v_end   = 4.5;                              # end

v_nodes = 100;                              # velocity discretization, velocity-points

v = collect(range(v_start,v_end,v_nodes));  # collect(range(a,b,n)) <-> LinRange(x): works!
w = 9/v_nodes*ones(size(v));                # width (for numerical integration)

########## calculation initial elements of U ##########

xx = LinRange(x_start,x_end,k+1);
xm = 0.5*(xx[1:end-1] + xx[2:end]);

U       = zeros(k,3);
g_val   = zeros(k+1,v_nodes);

U[:,1] = (xm.<=0.5)*Prim_l[1]+(xm.>0.5)*Prim_r[1];      # rho
U[:,2] = (xm.<=0.5)*Prim_l[2]+(xm.>0.5)*Prim_r[2];      # u bulk velocity of the gas
U[:,3] = (xm.<=0.5)*Prim_l[3]+(xm.>0.5)*Prim_r[3];      # temperatur

########## initial data: Flux, Source ##########

Flux    = zeros(k+1,3);
Source  = zeros(k,3);

########## parts of v ##########

vvp     = kron(ones(k+1,1),v'.*(v'.>0));    # v^{+}
vvm     = kron(ones(k+1,1),v'.*(v'.<0));    # v^{-}
vv      = kron(ones(k+1,1),v');             # v
vv_s    = kron(ones(k,1),v');

dt      = 0.0005;                           # time-step
t_end   = 0.14;                             # running time

while t < t_end

M               = Maxwell(U, v);
vM              = M[[1;1:k],:].*vvp + M[[1:k;k],:].*vvm;                                                        # term: <m*(v^{+}*M^{n}_{i} + v^{-}*M^{n}_{i+1})>

vg_dx           = (g_val[1:k+1,:]-g_val[[1;1:k],:]).*vvp/dx + (g_val[[2:k+1;k+1],:]-g_val[1:k+1,:]).*vvm/dx;    # term: [v^{+}*(g^{n}_{i+1/2}-g^{n}_{i-1/2})/deltax + v^{-}*(g^{n}_{i+3/2}-g^{n}_{i+1/2})/deltax]
vM_dx           = vv.*(M[[1:k;k],:]-M[[1;1:k],:])/dx;                                                           # term: (v*(M^{n}_{i+1}-M^{n}_{i})/deltax)

A               = id_min_proj(U, M, vg_dx, vM_dx, v, w);                        # return value id_min_proj: tensor([251,100],2)
vg_dx           = A[1];
vM_dx           = A[2];

g_val           = (g_val - dt*vg_dx - dt/epsilon*vM_dx)/(1+dt/epsilon);         # kinetic equation

Flux[:,1]       = vM*w;
Flux[:,2]       = (vv.*vM)*w;
Flux[:,3]       = 0.5*(vv.^2 .*vM)*w;

g_dx            = (g_val[2:k+1,:]-g_val[1:k,:])/dx;
Source[:,1]     = g_dx*diagm(w)*v;
Source[:,2]     = (vv_s.*g_dx)*diagm(w)*v;
Source[:,3]     = 0.5*(vv_s.^2 .*g_dx)*diagm(w)*v;

U               = U - dt/dx*(Flux[2:end,:]-Flux[1:end-1,:]) - dt*epsilon*Source;        # macroscopic equation

t               = t+dt;

end

rho = U[:,1];
u   = U[:,2]./rho;
T   = 2*U[:,3]./rho-u.^2;
``````
``````using LinearAlgebra

function id_min_proj(U, M, vg, vM, v, w)

rho = [U[1,1];U[:,1]];
u   = [U[2,2];U[:,2]]./rho;
T   = 2*[U[3,3];U[:,3]]./rho-u.^2;

M = [M[1,:]';M];

vv = kron(ones(size(rho)),v');
uu = kron(u,ones(size(v')));
TT = kron(T,ones(size(v')));
rrho = kron(rho,ones(size(v')));

########## terms of vg ##########

mean_phi    = vg*w;                                     # term: <phi>
mean_vphi   = ((vv-uu).*vg)*w;                          # term: <(v-u)*phi>
mean_v2phi  = ((0.5*(vv-uu).^2 ./TT .-0.5).*vg)*w;      # term: <((v-u)^2/2T - 1/2)*phi>; d=1, because one-dimensional gas

vg_proj = (kron((mean_phi./rho), ones(size(v'))) + kron((mean_vphi./(rho.*T)), ones(size(v'))).*(vv-uu) + kron((mean_v2phi./rho), ones(size(v'))).*((vv-uu).^2 ./TT .-1)).*M;
ret_vg = vg_proj;

########## terms of vM ##########

mean_phi    = vM*w;                                     # term: <phi>
mean_vphi   = ((vv-uu).*vM)*w;                          # term: <(v-u)*phi>
mean_v2phi  = ((0.5*(vv-uu).^2 ./TT .-0.5).*vM)*w;      # term: <((v-u)^2/2T - 1/2)*phi>; d=1, because one-dimensional gas

vM_proj = (kron((mean_phi./rho), ones(size(v'))) + kron((mean_vphi./(rho.*T)), ones(size(v'))).*(vv-uu)  + kron((mean_v2phi./rho), ones(size(v'))).*((vv-uu).^2 ./TT .-1)).*M;
ret_vM = vM_proj;

return ret_vg, ret_vM

end
``````
``````using LinearAlgebra

function Maxwell(U, v)

rho     = U[:,1];                       # rho
u       = U[:,2]./rho;                  # u bulk velocity of the gas
T       = 2*U[:,3]./rho-u.^2;           # temperatur

vv      = kron(ones(size(rho)),v');
uu      = kron(u,ones(size(v')));
TT      = kron(T,ones(size(v')));
rrho    = kron(rho,ones(size(v')));

return M = rrho./sqrt.(2*pi*TT).*exp.((-0.5*(vv-uu).^2)./TT);

end

``````

I get the following output in REPL:

``````┌ Warning: Assignment to `g_val` in soft scope is ambiguous because a global variable by the same name exists: `g_val` will be treated as a new local. Disambiguate by using `local g_val` to suppress this warning or `global g_val` to assign to the existing global variable.
└ @ ~/Julia.Projekte/Test.jl:71
┌ Warning: Assignment to `U` in soft scope is ambiguous because a global variable by the same name exists: `U` will be treated as a new local. Disambiguate by using `local U` to suppress this warning or `global U` to assign to the existing global variable.
└ @ ~/Julia.Projekte/Test.jl:82
┌ Warning: Assignment to `t` in soft scope is ambiguous because a global variable by the same name exists: `t` will be treated as a new local. Disambiguate by using `local t` to suppress this warning or `global t` to assign to the existing global variable.
└ @ ~/Julia.Projekte/Test.jl:84
ERROR: UndefVarError: `U` not defined
Stacktrace:
[1] top-level scope
@ ~/Julia.Projekte/Test.jl:61

┌ Warning: Assignment to `g_val` in soft scope is ambiguous because a global variable by the same name exists: `g_val` will be treated as a new local. Disambiguate by using `local g_val` to suppress this warning or `global g_val` to assign to the existing global variable.
└ @ ~/Julia.Projekte/Test.jl:71
┌ Warning: Assignment to `U` in soft scope is ambiguous because a global variable by the same name exists: `U` will be treated as a new local. Disambiguate by using `local U` to suppress this warning or `global U` to assign to the existing global variable.
└ @ ~/Julia.Projekte/Test.jl:82
┌ Warning: Assignment to `t` in soft scope is ambiguous because a global variable by the same name exists: `t` will be treated as a new local. Disambiguate by using `local t` to suppress this warning or `global t` to assign to the existing global variable.
└ @ ~/Julia.Projekte/Test.jl:84
ERROR: UndefVarError: `U` not defined
Stacktrace:
[1] top-level scope
@ ~/Julia.Projekte/Test.jl:61
``````

Okay … the matrix M is not calculated, I think, due to the warnings. I took some looks at books about Julia with code examples. I cannot find any hints with this scopes. All in all … this code works in Matlab. I have re-write this code in Julia, with some changes (such as the matrix-notations).

The code seems to be work, when I set’ global’ (but I think, this is the the final BEST solution!):

``````using LinearAlgebra
using Plots

t = 0;

epsilon = 0.1;

########## initial data: Sod problem (rho,u,T) ##########

Prim_l = [1 0 1];                           # 0 <= x <= 0.5
Prim_r = [0.125 0 0.1];                     # 0.5 < x <= 1

########## domain in space ##########

x_start = 0;                                # start
x_end   = 1;                                # end

k   = 250;                                  # cell discretization, grid-points
dx  = (x_end-x_start)/k;

########## velocity space (v = molecular velocities at certain point x at time t) ##########

v_start = -4.5;                             # start
v_end   = 4.5;                              # end

v_nodes = 100;                              # velocity discretization, velocity-points

v = collect(range(v_start,v_end,v_nodes));  # collect(range(a,b,n)) <-> LinRange(x): works!
w = 9/v_nodes*ones(size(v));                # width (for numerical integration)

########## calculation initial elements of U ##########

xx = LinRange(x_start,x_end,k+1);
xm = 0.5*(xx[1:end-1] + xx[2:end]);

U       = zeros(k,3);
g_val   = zeros(k+1,v_nodes);

U[:,1] = (xm.<=0.5)*Prim_l[1]+(xm.>0.5)*Prim_r[1];      # rho
U[:,2] = (xm.<=0.5)*Prim_l[2]+(xm.>0.5)*Prim_r[2];      # u bulk velocity of the gas
U[:,3] = (xm.<=0.5)*Prim_l[3]+(xm.>0.5)*Prim_r[3];      # temperatur

########## initial data: Flux, Source ##########

Flux    = zeros(k+1,3);
Source  = zeros(k,3);

########## parts of v ##########

vvp     = kron(ones(k+1,1),v'.*(v'.>0));    # v^{+}
vvm     = kron(ones(k+1,1),v'.*(v'.<0));    # v^{-}
vv      = kron(ones(k+1,1),v');             # v
vv_s    = kron(ones(k,1),v');

dt      = 0.0005;                           # time-step
t_end   = 0.14;                             # running time

while t < t_end

global M               = Maxwell(U, v);
global vM              = M[[1;1:k],:].*vvp + M[[1:k;k],:].*vvm;                                                        # term: <m*(v^{+}*M^{n}_{i} + v^{-}*M^{n}_{i+1})>

global vg_dx           = (g_val[1:k+1,:]-g_val[[1;1:k],:]).*vvp/dx + (g_val[[2:k+1;k+1],:]-g_val[1:k+1,:]).*vvm/dx;    # term: [v^{+}*(g^{n}_{i+1/2}-g^{n}_{i-1/2})/deltax + v^{-}*(g^{n}_{i+3/2}-g^{n}_{i+1/2})/deltax]
global vM_dx           = vv.*(M[[1:k;k],:]-M[[1;1:k],:])/dx;                                                           # term: (v*(M^{n}_{i+1}-M^{n}_{i})/deltax)

global A               = id_min_proj(U, M, vg_dx, vM_dx, v, w);                        # return value id_min_proj: tensor([251,100],2)
vg_dx           = A[1];
vM_dx           = A[2];

global g_val           = (g_val - dt*vg_dx - dt/epsilon*vM_dx)/(1+dt/epsilon);         # kinetic equation

Flux[:,1]       = vM*w;
Flux[:,2]       = (vv.*vM)*w;
Flux[:,3]       = 0.5*(vv.^2 .*vM)*w;

global g_dx            = (g_val[2:k+1,:]-g_val[1:k,:])/dx;
Source[:,1]     = g_dx*diagm(w)*v;
Source[:,2]     = (vv_s.*g_dx)*diagm(w)*v;
Source[:,3]     = 0.5*(vv_s.^2 .*g_dx)*diagm(w)*v;

global U               = U - dt/dx*(Flux[2:end,:]-Flux[1:end-1,:]) - dt*epsilon*Source;        # macroscopic equation

global t               = t+dt;

end

rho = U[:,1];
u   = U[:,2]./rho;
T   = 2*U[:,3]./rho-u.^2;

plot(xm, rho)
``````

Kind regars!

Hi there!
You kinds disregarded the most important rule in Julia: don’t compute stuff in global scope! It’s slow and the scoping works a bit differently (as you experience).

So the quick fix is just to wrap your whole code (expect the `using` statements) in the first code block you posted in a function (e.g. `main()`) and then call it. That should be sufficient to make everything work (if the code is otherwise correct).

Hello,

When were scopes introduced in Julia?
I have some books about Julia, and nothing is mentioned about scopes there.
I used some code example from ‘Ferreira Introduction to Computational Physics with examples in Julia, 2016’.
Of course, these codes work with Julia due to scope problems.

I can’t tell you when the current scoping rules where implemented, but I think the latest change was soft global scope for the REPL which should’nt influence your code. You can read about scopes here:
https://docs.julialang.org/en/v1/manual/variables-and-scoping/

Please note that examples from 2016 are ancient with respect to Julia. Julia 1.0 was released in 2018. So you should probably find a more upto date source.

How does the Main() function look like, for example? I need corresponding inputs for functions, which are defined in brackets.

The quick fix is really just wrapping everything (except `using` statements) in a function and then calling that:

``````using LinearAlgebra
using Plots

function main() # wrap everything in a function
t = 0;

epsilon = 0.1;

########## initial data: Sod problem (rho,u,T) ##########

Prim_l = [1 0 1];                           # 0 <= x <= 0.5
Prim_r = [0.125 0 0.1];                     # 0.5 < x <= 1

########## domain in space ##########

x_start = 0;                                # start
x_end   = 1;                                # end

k   = 250;                                  # cell discretization, grid-points
dx  = (x_end-x_start)/k;

########## velocity space (v = molecular velocities at certain point x at time t) ##########

v_start = -4.5;                             # start
v_end   = 4.5;                              # end

v_nodes = 100;                              # velocity discretization, velocity-points

v = collect(range(v_start,v_end,v_nodes));  # collect(range(a,b,n)) <-> LinRange(x): works!
w = 9/v_nodes*ones(size(v));                # width (for numerical integration)

########## calculation initial elements of U ##########

xx = LinRange(x_start,x_end,k+1);
xm = 0.5*(xx[1:end-1] + xx[2:end]);

U       = zeros(k,3);
g_val   = zeros(k+1,v_nodes);

U[:,1] = (xm.<=0.5)*Prim_l[1]+(xm.>0.5)*Prim_r[1];      # rho
U[:,2] = (xm.<=0.5)*Prim_l[2]+(xm.>0.5)*Prim_r[2];      # u bulk velocity of the gas
U[:,3] = (xm.<=0.5)*Prim_l[3]+(xm.>0.5)*Prim_r[3];      # temperatur

########## initial data: Flux, Source ##########

Flux    = zeros(k+1,3);
Source  = zeros(k,3);

########## parts of v ##########

vvp     = kron(ones(k+1,1),v'.*(v'.>0));    # v^{+}
vvm     = kron(ones(k+1,1),v'.*(v'.<0));    # v^{-}
vv      = kron(ones(k+1,1),v');             # v
vv_s    = kron(ones(k,1),v');

dt      = 0.0005;                           # time-step
t_end   = 0.14;                             # running time

while t < t_end

M               = Maxwell(U, v);
vM              = M[[1;1:k],:].*vvp + M[[1:k;k],:].*vvm;                                                        # term: <m*(v^{+}*M^{n}_{i} + v^{-}*M^{n}_{i+1})>

vg_dx           = (g_val[1:k+1,:]-g_val[[1;1:k],:]).*vvp/dx + (g_val[[2:k+1;k+1],:]-g_val[1:k+1,:]).*vvm/dx;    # term: [v^{+}*(g^{n}_{i+1/2}-g^{n}_{i-1/2})/deltax + v^{-}*(g^{n}_{i+3/2}-g^{n}_{i+1/2})/deltax]
vM_dx           = vv.*(M[[1:k;k],:]-M[[1;1:k],:])/dx;                                                           # term: (v*(M^{n}_{i+1}-M^{n}_{i})/deltax)

A               = id_min_proj(U, M, vg_dx, vM_dx, v, w);                        # return value id_min_proj: tensor([251,100],2)
vg_dx           = A[1];
vM_dx           = A[2];

g_val           = (g_val - dt*vg_dx - dt/epsilon*vM_dx)/(1+dt/epsilon);         # kinetic equation

Flux[:,1]       = vM*w;
Flux[:,2]       = (vv.*vM)*w;
Flux[:,3]       = 0.5*(vv.^2 .*vM)*w;

g_dx            = (g_val[2:k+1,:]-g_val[1:k,:])/dx;
Source[:,1]     = g_dx*diagm(w)*v;
Source[:,2]     = (vv_s.*g_dx)*diagm(w)*v;
Source[:,3]     = 0.5*(vv_s.^2 .*g_dx)*diagm(w)*v;

U               = U - dt/dx*(Flux[2:end,:]-Flux[1:end-1,:]) - dt*epsilon*Source;        # macroscopic equation

t               = t+dt;

end

rho = U[:,1];
u   = U[:,2]./rho;
T   = 2*U[:,3]./rho-u.^2;
end

main() # now call it
``````