Domains and Material Density#

Real micromagnetic problems are rarely a single uniform block. You usually want a disk in vacuum, a stack of two materials, a multi-layer with an air spacer, or a topology-optimized geometry where the “material vs. air” question is itself the design variable. NeuralMag handles all of this through a single concept — the material density state.rho — and a small set of helpers built on top of it.

What rho is#

state.rho is a CellFunction that takes the value 1.0 inside the magnetic material and state.eps (machine epsilon) outside. It is created automatically when you instantiate a State:

# neuralmag/common/state.py
self.rho = CellFunction(
    self, tensor=lambda domain: np.where(domain, 1.0, self.eps)
)

So rho is a dynamic attribute of the state (see Dynamic Attributes, resolve and remap): it depends on the boolean field state.domain, which in turn depends on state.domains (the integer domain-id field). You will normally not touch the lambda — you only populate state.domains via State.add_domain() and let the rest follow.

Why eps and not zero?

Several field terms divide by mass-lumped denominators that contain rho. A density of exactly zero would produce inf / nan instead of “no contribution from this cell”. Use state.eps (or any other small positive number) for non-magnetic regions.

How rho enters the energy#

Every effective-field and energy form generated by the form compiler multiplies its integrand by rho. This is why you can write energy expressions like

return A * dot(grad(m), grad(m)) * dV(dim)

without any explicit mask: cells with rho ≈ eps simply drop out of the integral. The same trick is what makes topology optimization differentiable — gradients with respect to rho flow through every field term automatically.

Interface densities (3D only)#

In 3D, the state additionally exposes three nodal density fields,

state.rhoxy   # ("ccn") density at xy-interfaces
state.rhoxz   # ("cnc") density at xz-interfaces
state.rhoyz   # ("ncc") density at yz-interfaces

These are generated automatically (see State._generate_code) by taking the maximum of the cell rho values on either side of each interface. Surface-energy terms (e.g. interface DMI) integrate against them so that an interface between a magnetic cell and an air cell still contributes, while two-air interfaces vanish. You normally don’t touch them by hand.

Defining regions with add_domain#

State.add_domain() assigns an integer ID to every cell satisfying a boolean condition. Domain 0 is the implicit background; the first call to add_domain zeroes out the default fill of 1 and starts marking cells.

# demos/skyrmion-disk.py
x, y, z = state.coordinates()
state.add_domain(1, x ** 2 + y ** 2 < 50e-9 ** 2)

After this call state.domains has value 1 inside the disk and 0 outside; state.rho automatically becomes 1.0 inside and state.eps outside.

You can layer several domains:

x, _, _ = state.coordinates()
L0 = 30e-9
state.add_domain(1, abs(x) <= L0)   # chiral region
state.add_domain(2, abs(x) >  L0)   # ferromagnetic slabs

Re-adding a domain ID overwrites previously marked cells in those locations, so the order of calls matters when conditions overlap.

Per-domain material parameters with fill_by_domain#

The most common reason to use multiple domains is to assign different material parameters per region. CellFunction.fill_by_domain() does exactly this: it takes a list of values, one per domain ID (starting at 0), and stores a tiny lambda that indexes into that list using state.domains.

import neuralmag as nm

x, _, _ = state.coordinates()
state.add_domain(1, abs(x) <= L0)   # chiral
state.add_domain(2, abs(x) >  L0)   # ferro

# values are indexed as [domain_0_background, domain_1_chiral, domain_2_ferro]
state.material.Ms = nm.CellFunction(state).fill_by_domain([0.0, Ms, Ms])
state.material.A  = nm.CellFunction(state).fill_by_domain([0.0, A,  A ])
state.material.Db = nm.CellFunction(state).fill_by_domain([0.0, D,  0.0])
state.material.Ku = nm.CellFunction(state).fill_by_domain([0.0, Kc, Ku])

Because fill_by_domain stores a lambda on the function, the values re-index automatically if you ever update state.domains. The function is currently restricted to fully cell-based spaces ("c" * dim); attempts to call it on a nodal function raise NotImplementedError.

Topology-optimization pattern#

For a true topology-optimization problem you do not enumerate domains — you make rho itself a function of an underlying continuous design variable rho_m:

# demos/topology-optimization_jax.py
state.rho   = nm.CellFunction(
    state,
    tensor=lambda rho_m: rho_m.at[:, :, 10:].set(state.eps),
)
state.rho_m = nm.CellFunction(state).fill(1.0)

Here state.rho is overridden to (a) take the design variable rho_m as input and (b) force the upper layers to be air. From this point on, state.resolve("h_demag", ["rho_m"]) (see Dynamic Attributes, resolve and remap) gives you a pure function of rho_m that you can plug into jax.grad.

Gotchas

  • Never set any cell’s rho to exactly zero — use state.eps.

  • fill_by_domain only works on fully cell-based functions.

  • Domain IDs must be non-negative integers; ID 0 is reserved for the background.

  • Re-adding the same ID overwrites earlier marks; if conditions overlap, the last call wins.

  • In 1D / 2D the interface densities rhoxy / rhoxz / rhoyz do not exist — they are 3D-only by construction.

See also#