How to define input from users for a 2D heteregenous model

Hi everyone,

I am working on a small modeling project. Let’s say I want to do 2D diffusion modeling in different units. All the units are homogeneous but they can have different geometries and physical properties.

I was wondering what would be the most Julia way to ask the user for that?

My first idea is to ask the user to give a raster image with each pixel color beeing a different units or to draw on a GUI his model. He would then need to give names to each units and tell me their physical properties.

From that, I can build boolean matrices for the position of each units and define matrices for each physical properties for the whole model.

What would be your idea for this situation?

The most flexible thing is probably to have the caller pass in functions f(x,y) rather than images. (This is called a higher-order function.)

e.g. suppose you are solving the diffusion equation du/dt = \nabla \cdot (D \nabla u) + f with some initial conditions u(x,y,0) = u_0(x,y) and some boundary conditions on a box domain [0,L_x] \times [0,L_y], say with finite differences. Then you ask the user to specify the size L_x \times L_y of the domain, and pass in functions f(x,y) and u_0(x,y) and D(x,y), along with some information about the desired discretization and the simulation time.

Note that your program doesn’t need to care about the units of L and f etcetera, as long as the caller was consistent in their units.


Mmh, I am a bit confused. Let’s say the user wants to have a geometry like this for his model:


with the three colors having different values of D, f and u0

How can he build that easily as a function?

The other thing is that my project is not targeting people who knows programming. So I can’t expect them to write any complicated functions.


Alright, after thinking a bit, I guess your idea would be for the user, for example, to write something like that for D(x,y)

function D_calc(nx,ny)

    D = zeros(nx,ny)

    D[1:Int(ny/10), :] .= 1
    D[Int(ny/10)+1:end, :] .= 2

    return D

to obtain:

and to use that function as an input for my model (correct me if I am wrong). But that is still a bit to complicated for the people I am targetting… And it should be a bit more tricky for more fancy geometries.

There’s really no good way to initialize simulations without any programming knowledge unless you’re willing to build a GUI, which can be quite expensive and time-consuming. Higher-order functions may work well enough if you provide example scripts and a few geometric primitives (box, circle, line, ellipsoid, …). I’d recommend Unitful.jl to handle different input units on the frontend, with everything converted internally to a consistent unit system on the backend.

As for specification of higher-order functions, the user should provide a function f(x, y) that’s independent of the grid. If all material interfaces are initially sharp (i.e. materials start unmixed), you could do something like this:

using Unitful

@derived_dimension MassDiffusivity  Unitful.𝐋^2/Unitful.𝐓

struct DiffusiveMaterial
    D::Float64  # Mass diffusivity
    ρ::Float64  # Density
    function DiffusiveMaterial(D, ρ)
        D < 0 && error("Diffusivity must be positive")
        ρ < 0 && error("Density must be positive")
        new(D, ρ)

# Convert to uniform internal units
function DiffusiveMaterial(D::MassDiffusivity, ρ::Unitful.Density)
    DiffusiveMaterial(ustrip(u"m^2/s", D), ustrip(u"kg/m^3", ρ))

struct Domain{T <: AbstractRange}
    D::Array{Float64, 2}
    ρ::Array{Float64, 2}
    materials::Dict{Symbol, DiffusiveMaterial}
    function Domain(x::T, y::T, materials) where {T}
        return new{T}(x, y, 
                      Array{Float64}(undef, length(x), length(y)), 
                      Array{Float64}(undef, length(x), length(y)), 

function initialize!(material, d::Domain)
    for (j, y) in enumerate(d.y), (i, x) in enumerate(d.x)
        m = material(x, y)
        d.D[i, j] = d.materials[m].D
        d.ρ[i, j] = d.materials[m].ρ

incircle(x₀, y₀, R) = (x, y) -> hypot(x - x₀, y - y₀) <= R
inrectangle(x₀, y₀, x₁, y₁) = (x, y) -> x₀ <= x <= x₁ && y₀ <= y <= y₁

which would allow the user to initialize a model like this

materials = Dict(:Mud =>       DiffusiveMaterial(1.0e-7u"cm^2/s", 2.65e3u"kg/m^3"),
                 :Water =>     DiffusiveMaterial(1.0e-6u"cm^2/s", 1000u"kg/m^3"),
                 :Ethanol =>   DiffusiveMaterial(8.0e-6u"cm^2/s", 800u"kg/m^3"),)

d = Domain(LinRange(0, 100, 512), LinRange(0, 200, 1024), materials)

initialize!(d) do x, y
    if incircle(50, 50, 20)(x, y)
    elseif inrectangle(0, 0, 200, 30)(x, y)

Thank you very much for taking the time to write that! I’ve learned a lot about structures and how I could make use of that by playing with your code. Very elegant!