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

General Usage
,
1

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?

2

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.

2 Likes
3

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

test

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.

EDIT

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
end

to obtain:
plot_2

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.

4

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, ρ)
    end
end

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

struct Domain{T <: AbstractRange}
    x::T
    y::T
    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)), 
                      materials)
    end
end

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].ρ
    end
end

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)
        :Ethanol
    elseif inrectangle(0, 0, 200, 30)(x, y)
        :Mud
    else
        :Water
    end
end
2 Likes
5

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!

2 Likes