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