Numerical-Code: Warnings and Errors due to scopes

Specific Domains Modelling & Simulations
1

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!

2

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

3

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.

4

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.

5

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

6

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