You can just define the function with a default parameter
func1(x; amp=DEFAULT_VALUE)
and pass a different value to it when needed
2 Likes
Thanks for the answer!
I want to develop a physics simulation package. Every physical quantity is normalized by a user-chosen characteristic length, says lambda.
module M
func1(t; lambda) = do_something
func2(t; lambda) = do_something
func3(t; lambda) = do_something
func4(t; lambda) = do_something
end
func1, func2, func3, func4 are part of simulation procedure. Since there will be only one characteristic length in the entire simulation, I am trying to provide a way to set the characteristic length before the simulation, instead of passing it every time when calling the simulation procedure. Would it be possible to achieve that?
Thanks for the answer!
Not sure if this approach has the similar performance issue mentioned by @stevengj ?
Make a Simulation struct that contains lambda and have functions operate on that object.
6 Likes
The performance would be similar as long as you get rid of global variables.
1 Like
I created the Jupyter and Pluto notebooks with the same contents as below:
- Jupyter notebooks:
- Pluto notebooks:
Very Long!
Problem-Algorithm-Solver pattern
A code will be more readable and efficient if the global variables that store the parameters of a problem are passed as arguments to functions everytime.
However, many people find it troublesome to pass the parameters as arguments to the function everytime.
In order to resolve this misunderstanding, they can try the problem-algorithm-solver pattern explained below.
In order to use the (; a, b, c) = p syntax, require VERSION ≥ v"1.7.0-beta" #39285.
VERSION ≥ v"1.7.0-beta" #> true
A minimal working example of the problem-algorithm-solver pattern:
module FreeFall
"""Problem type"""
Base.@kwdef struct Problem{G, Y0, V0, TS}
g::G = 9.80665
y0::Y0 = 0.0
v0::V0 = 30.0
tspan::TS = (0.0, 8.0)
end
"""Algorithm types"""
Base.@kwdef struct EulerMethod{T} dt::T = 0.1 end
Base.@kwdef struct ExactFormula{T} dt::T = 0.1 end
default_algorithm(::Problem) = EulerMethod()
"""Solution type"""
struct Solution{Y, V, T, P<:Problem, A} y::Y; v::V; t::T; prob::P; alg::A end
"""Solver function"""
solve(prob::Problem) = solve(prob, default_algorithm(prob))
function solve(prob::Problem, alg::ExactFormula)
(; g, y0, v0, tspan) = prob
(; dt) = alg
t0, t1 = tspan
t = range(t0, t1 + dt/2; step = dt)
y(t) = y0 + v0*(t - t0) - g*(t - t0)^2/2
v(t) = v0 - g*(t - t0)
Solution(y.(t), v.(t), t, prob, alg)
end
function solve(prob::Problem, alg::EulerMethod)
(; g, y0, v0, tspan) = prob
(; dt) = alg
t0, t1 = tspan
t = range(t0, t1 + dt/2; step = dt)
n = length(t)
y = Vector{typeof(y0)}(undef, n)
v = Vector{typeof(v0)}(undef, n)
y[1] = y0
v[1] = v0
for i in 1:n-1
v[i+1] = v[i] - g*dt
y[i+1] = y[i] + v[i]*dt
end
Solution(y, v, t, prob, alg)
end
end
This module will calculate the free-fall motion under a uniform gravity potential.
-
Problem: The parameters of problems are stored in objects of type
FreeFall.Problem. -
Algorithm: The parameters of algorithms are stored in objects of types
FreeFall.EulerMethodand ofFreeFall.ExactFormula. -
Solver function: The function
FreeFall.solve(prob::FreeFall.Problem, alg)returns the object of typeFreeFall.Solutionobtained by solving the problemprobwith the algorithmalg.
Note that the (; a, b, c) = p syntax of Julia ≥ v1.7, is very useful for extracting the contents of the object that contains parameters. (You can also do the same thing with @unpack in the exellent well-known package Parameters.jl.)
Example 1: Compare a numerical solution and an exact one.
using Plots
earth = FreeFall.Problem()
sol_euler = FreeFall.solve(earth)
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())
plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)
Since the time step dt = 0.1 is not enough small, the error of the Euler method is rather large.
Example 2: Change the algorithm parameter dt smaller.
earth = FreeFall.Problem()
sol_euler = FreeFall.solve(earth, FreeFall.EulerMethod(dt = 0.01))
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())
plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)
Example 3: Change the problem parameter g to the gravitational acceleration on the moon.
moon = FreeFall.Problem(g = 1.62, tspan = (0.0, 40.0))
sol_euler = FreeFall.solve(moon)
sol_exact = FreeFall.solve(moon, FreeFall.ExactFormula(dt = sol_euler.alg.dt))
plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Moon"; xlabel="t", legend=:bottomleft)
Extensibility of the Julia ecosystem: A different author of the FreeFall module can define , in the SomeExtension module, a new algorithm type and a new method of the FreeFall.solve function. Both a new type and a new method!
module SomeExtension
using ..FreeFall
using ..FreeFall: Problem, Solution
Base.@kwdef struct Symplectic2ndOrder{T} dt::T = 0.1 end
function FreeFall.solve(prob::Problem, alg::Symplectic2ndOrder)
(; g, y0, v0, tspan) = prob
(; dt) = alg
t0, t1 = tspan
t = range(t0, t1 + dt/2; step = dt)
n = length(t)
y = Vector{typeof(y0)}(undef, n)
v = Vector{typeof(v0)}(undef, n)
y[1] = y0
v[1] = v0
for i in 1:n-1
ytmp = y[i] + v[i]*dt/2
v[i+1] = v[i] - g*dt
y[i+1] = ytmp + v[i+1]*dt/2
end
Solution(y, v, t, prob, alg)
end
end
earth = FreeFall.Problem()
sol_sympl = FreeFall.solve(earth, SomeExtension.Symplectic2ndOrder(dt = 2.0))
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())
plot(sol_sympl.t, sol_sympl.y; label="2nd order symplectic (dt = $(sol_sympl.alg.dt))", ls=:auto)
plot!(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft)
The second-order symplectic method gives exact discretizations of the free-fall motions under uniform gravity potentials.
Composability of the Julia ecosystem: By combining the MonteCarloMeasurements.jl package with the modules given above, we can plot the behavior of the system for perturbations of the initial value. We did not intend for them to be used in this way!
using MonteCarloMeasurements
earth = FreeFall.Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)
sol_euler = FreeFall.solve(earth)
sol_sympl = FreeFall.solve(earth, SomeExtension.Symplectic2ndOrder(dt = 2.0))
sol_exact = FreeFall.solve(earth, FreeFall.ExactFormula())
ylim = (-100, 60)
P = plot(sol_euler.t, sol_euler.y; label="Euler's method (dt = $(sol_euler.alg.dt))", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
Q = plot(sol_sympl.t, sol_sympl.y; label="2nd order symplectic (dt = $(sol_sympl.alg.dt))", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
R = plot(sol_exact.t, sol_exact.y; label="exact solution", ls=:auto)
title!("On the Earth"; xlabel="t", legend=:bottomleft, ylim)
plot(P, Q, R; size=(720, 600))
Edit 1: Typo. PLots → Plots
Edit 2: Add Jupyter notebook
Edit 3: Add Pluto notebooks and Julia v1.6 versions
15 Likes
Globals make your code unusable in multithreading. If someone tries to solve a physics problem with your package in parallel for different lengths and you use globals… It will Bork heavily. Use the suggested approaches above and it will be possible to multithread!
4 Likes
Thanks for your answer. Although my case is easier than the example you provided, I learn a lot from your example.
What is the purpose of using the 4 different type parameters (G, Y0, V0, TS) in the declaration of Problem structure?
struct Problem{G, Y0, V0, TS}
g::G = 9.80665
y0::Y0 = 0.0
v0::V0 = 30.0
tspan::TS = (0.0, 8.0)
end
Why not just using one type parameter instead
struct Problem{T}
g::T = 9.80665
y0::T = 0.0
v0::T = 30.0
tspan::Tuple{T, T} = (0.0, 8.0)
end
1 Like
I created the Jupyter notebook with this change:
Most results are the same as the original. However, only the last example results in an error:
MethodError: no method matching Main.FreeFall.Problem(::Float64, ::Particles{Float64, 2000}, ::Particles{Float64, 2000}, ::Tuple{Float64, Float64})
Closest candidates are:
Main.FreeFall.Problem(::T, ::T, ::T, ::Tuple{T, T}) where T at In[2]:13
In the last example of using MonteCarloMeasurements, the types of the fields y0 and v0 are set to MonteCarloMeasurements.Particles{Float64, 2000} instead of Float64.
Fixed version:
Base.@kwdef struct Problem{T, U}
g::T = 9.80665
y0::U = 0.0
v0::U = 30.0
tspan::Tuple{T, T} = (0.0, 8.0)
end
Thus, overly restricting the type combinations may inhibit the composability, the strengths of the Julia ecosystem.
When in doubt, it is safe to mechanically separate the type parameters of all potentially different types of fields.
Doing so has the advantage of reducing the amount of something to worry about and makes it easier to concentrate on mathematics and physics. 
2 Likes
@genkuroki, thanks for this outstanding post, a great service to promote Julia language IMHO.
As an exercise, added to your code the different physical units using Unitful and UnitfulRecipes. It worked everywhere at the first attempt except for the MonteCarloMeasurements plots. Do you know of any recipe for that?
NB: will be glad to post the preliminary edited code if useful.
By saying multithreading, do you mean things like openMP?
However, if use MPI instead, will global variables still have such problems you mentioned?
Is it always recommended to parameterize the type of a struct? For instance, if a struct has 10 fields, then do we need to specify 10 type parameters for each field?
Sorry for a stupid question, what does the … mean in
using ..FreeFall
I know if I write a module called FreeFall.jl, I can do
include("FreeFall.jl")
using .FreeFall
But what does the … do in using …FreeFall?
Thanks!
Yes.
First, for type stability. (cf. Avoid fields with abstract type)
Second, for composability with various packages that cannot be predicted in advance. (Example: Composability with MonteCarloMeasurements.jl)
If you really don’t want to write type parameters, you can use ConcreteStructs.jl. The macro ConcreteStructs.@concrete automatically add type parameters to all fields in struct.
See [ANN] ConcreteStructs.jl - Cut the boilerplate when concretely parameterizing structs
Small working example. (Note that we have an issue someone should solve. Play nicely with other struct def macros · Issue #4 · jonniedie/ConcreteStructs.jl · GitHub)
using ConcreteStructs
Base.@kwdef @concrete struct Problem g; y0; v0; tspan end
Problem(;
g = 9.80665,
y0 = 0.0,
v0 = 30.0,
tspan = (0.0, 8.0)
) = Problem(g, y0, v0, tspan)
prob1 = Problem()
Output: Problem{Float64, Float64, Float64, Tuple{Float64, Float64}}(9.80665, 0.0, 30.0, (0.0, 8.0))
prob2 = Problem(v0 = 100.0)
Output: Problem{Float64, Float64, Float64, Tuple{Float64, Float64}}(9.80665, 0.0, 100.0, (0.0, 8.0))
using MonteCarloMeasurements
prob3 = Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)
Output: Problem{Float64, Particles{Float64, 2000}, Particles{Float64, 2000}, Tuple{Float64, Float64}}(9.80665, 0.0, 30.0 ± 1.0, (0.0, 8.0))
1 Like
Often you don’t want to, and instead want to parametrize many of your vars with the same parameter. For example,
Base.@kwdef struct Problem{T}
g::T = 9.80665
y0::T = 0.0
v0::T= 30.0
tspan::Tuple{T,T} = (0.0, 8.0)
end
2 Likes
Also, only if you need them to be parametric. They can also have static types:
struct A
x::Int
y::Float64
end
in which case you don’t need any parameter whatsoever.
2 Likes
no I mean things like
@threads for i in 1:n .... or @spawn begin .... end
Yes, that is basically “like OpenMP:” shared-memory threads, with decorator directives to make it easy. As opposed to distributed-memory programming like MPI.
Thanks. The second point above seems to be a design choice. Sometimes I feel that it is quite cumbersome to specify the argument type in a function if that type itself is a parametric one.
Consider again your above example. In solve(prob::Problem, alg::EulerMethod), no constraints are exerted on Problem. However, I usually have some constraints in mind during design. For example, I will not expect that the type parameter G to be a String (at least in my design). Then, shall I enforce this constraint in the signature of solve(prob::Problem, alg::EulerMethod) or simply throw an error during execution?
Thank you!
Yes, I was thinking it must be something related with shared memory stuff that can cause issue for global variable. For MPI, it seems using global variable is safe because each core are independent.