# Numerical-Code: Warnings and Errors due to scopes

**URL:** https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946
**Category:** Modelling & Simulations
**Created:** [August 8, 2024, 8:04am UTC](https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946 "2024-08-08T08:04:22Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![EatMyShorts](https://avatars.discourse-cdn.com/v4/letter/e/e5b9ba/32.png) [@EatMyShorts](https://discourse.julialang.org/u/EatMyShorts)
#### Post date: [August 8, 2024, 8:04am UTC](https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946/1 "2024-08-08T08:04:22Z")

</div>

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

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

```

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

```

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

```julia
┌ 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!):

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

---

<div class="post-metadata">

### Author: ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [August 10, 2024, 6:36am UTC](https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946/2 "2024-08-10T06:36:41Z")

</div>

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).

---

<div class="post-metadata">

### Author: ![EatMyShorts](https://avatars.discourse-cdn.com/v4/letter/e/e5b9ba/32.png) [@EatMyShorts](https://discourse.julialang.org/u/EatMyShorts)
#### Post date: [August 10, 2024, 10:22am UTC](https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946/3 "2024-08-10T10:22:07Z")

</div>

Hello,

thank you for your feedback!  
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.

---

<div class="post-metadata">

### Author: ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [August 10, 2024, 10:33am UTC](https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946/4 "2024-08-10T10:33:37Z")

</div>

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/](https://docs.julialang.org/en/v1/manual/variables-and-scoping/)

> [@EatMyShorts](#):
>
> with examples in Julia, 2016

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.

---

<div class="post-metadata">

### Author: ![EatMyShorts](https://avatars.discourse-cdn.com/v4/letter/e/e5b9ba/32.png) [@EatMyShorts](https://discourse.julialang.org/u/EatMyShorts)
#### Post date: [August 10, 2024, 7:38pm UTC](https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946/5 "2024-08-10T19:38:35Z")

</div>

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

---

<div class="post-metadata">

### Author: ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [August 10, 2024, 8:54pm UTC](https://discourse.julialang.org/t/numerical-code-warnings-and-errors-due-to-scopes/117946/6 "2024-08-10T20:54:08Z")

</div>

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

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

```