From Echoes to MeanFieldHomogenization

MeanFieldHomogenization (MFH) is a Julia port of Echoes [73], the C++ mean-field homogenization library (with a Python interface) developed at Cerema. If you already have Echoes scripts, this page is a direct translation guide: the same RVE/scheme/homogenize concepts, a line-by-line API table, and a worked example that calls both implementations — Echoes live from Julia via PyCall.jl — to check that they agree.

The @example blocks run at build time; the Python and PyCall blocks are shown for reference, and the Echoes numbers below were captured by running them against Echoes 1.0.

API translation table

Both libraries share the same mental model: an RVE holding phases, each with a geometry and a property dictionary, homogenized by a named scheme. They differ on one point besides syntax. Echoes names the matrix when the RVE is built, rve(matrix="SOLID"); MeanFieldHomogenization does not, because the same microstructure is a matrix/inclusion composite under Mori–Tanaka and a matrix-free aggregate under the self-consistent scheme. The phase that takes up the volume complement is declared fraction = :rest, and a matrix-based scheme picks it up unless it names another — so a translated script behaves the same.

Echoes (Python)	MeanFieldHomogenization (Julia)
`from echoes import *`	`using MeanFieldHomogenization, TensND`
`rve(matrix="SOLID")`	`RVE()`
`myrve["SOLID"] = ellipsoid(shape=spheroidal(1.), symmetrize=[ISO], prop={"C": stiff_kmu(k, μ)})`	`add_phase!(rve, :SOLID, Spheroid(1.0), Dict(:C => iso_stiffness(k, μ)); fraction = :rest, symmetrize = IsoSymmetrize())`
`myrve["PORE"] = ellipsoid(...)` then `myrve["PORE"].fraction = φ`	`add_phase!(rve, :PORE, Spheroid(1.0), Dict(:C => ...); fraction = φ, symmetrize = IsoSymmetrize())`
`crack(shape=spheroidal(ω), density=d, symmetrize=[ISO], prop={"C": tZ4})`	`add_phase!(rve, :CRACK, PennyCrack(1.0), Dict(:C => C0); density = d, symmetrize = IsoSymmetrize())`
`stiff_kmu(k, μ)`	`iso_stiffness(k, μ)`
`stiff_Enu(E, ν)`	`iso_stiffness_E_nu(E, ν)`
`tZ4` (zero 4th-order tensor, for voids/cracks)	`iso_stiffness(1e-6, 1e-6)` (near-zero, kept invertible)
`VOIGT`, `REUSS`	`Voigt`, `Reuss`
`DIL`, `DILD`	`Dilute`, `DiluteDual`
`MT`	`MoriTanaka`
`SC`, `ASC`	`SelfConsistent`, `AsymmetricSelfConsistent`
`MAX`, `PCW`	`Maxwell`, `PonteCastanedaWillis`
`DIFF`	`DifferentialScheme`
`homogenize(prop="C", rve=myrve, scheme=MT, epsrel=1e-6, maxnb=300, select_best=True)`	`homogenize(rve, MoriTanaka(), :C; reltol = 1e-6, abstol = 0, maxiters = 300, select_best = true)`
`C.k`, `C.mu`	`k_mu(C)`
`np.trace(K.array) / 3.`	`tr(Array(homogenize(rve, scheme, :K))) / 3`
`.paramsym(sym=TI)`	`best_fit_ti`

Two conventions worth flagging for a smooth transition:

Both stiff_kmu/iso_stiffness take physical (k, μ). Where they differ is in what the resulting tensor stores: MFH's raw TensISO{3}(a, b) constructor takes the pair (3k, 2μ), not (k, μ) — see the first tutorial. Building with iso_stiffness(k, μ) (as in the table above) sidesteps this entirely.
Symbol-string vs. type-instance schemes. Echoes selects a scheme via a module-level constant (MT, SC, …) passed to homogenize; MFH uses a scheme type instance (MoriTanaka(), SelfConsistent(; kwargs...)), which is also how solver options (abstol, reltol, maxiters, select_best) attach directly to the scheme rather than as loose homogenize keywords. A Symbol shortcut (:mt, :sc, …) is also accepted — see the schemes manual for the full alias table.

Same problem, both sides

The classic porous benchmark — a solid matrix with spherical pores, porosity — makes the translation concrete. In Echoes (the classic porous benchmark from the Echoes book):

python

from echoes import *

ks, mus = 72., 32.
kp, mup = 1.e-6, 1.e-6

myrve = rve(matrix="SOLID")
myrve["SOLID"] = ellipsoid(shape=spheroidal(1.), symmetrize=[ISO],
                            prop={"C": stiff_kmu(ks, mus)})
myrve["PORE"] = ellipsoid(shape=spheroidal(1.), symmetrize=[ISO],
                           prop={"C": stiff_kmu(kp, mup)})

def Chom_porous(phi, scheme):
    myrve["PORE"].fraction = phi
    myrve["SOLID"].fraction = 1. - phi
    C = homogenize(prop="C", rve=myrve, scheme=scheme,
                   epsrel=1.e-6, maxnb=300, select_best=True)
    return max(C.k, 0.), max(C.mu, 0.)

The same problem, live, in MFH:

julia

using MeanFieldHomogenization
using TensND
using LinearAlgebra
using Plots
gr()  # headless backend; GKSwstype is set to "100" in make.jl

const ks, mus = 72.0, 32.0
const kp, mup = 1.0e-6, 1.0e-6
C_s = iso_stiffness(ks, mus)
C_p = iso_stiffness(kp, mup)

function build_rve(φ)
    r = RVE()
    add_phase!(r, :SOLID, Spheroid(1.0), Dict(:C => C_s); fraction = :rest, symmetrize = IsoSymmetrize())
    add_phase!(r, :PORE, Spheroid(1.0), Dict(:C => C_p); fraction = φ, symmetrize = IsoSymmetrize())
    return r
end

φs = collect(0.0:0.1:1.0)
k_mt = [k_mu(homogenize(build_rve(φ), MoriTanaka(), :C))[1] for φ in φs]
k_mt[1:4]

4-element Vector{Float64}:
 72.0
 55.44385079613371
 43.06542136824175
 33.460581867908594

Line for line, the two snippets build the same RVE and call the same scheme — this is the translation to internalize for every other Echoes script.

Calling Echoes from Julia with PyCall

To cross-check MFH against a live Echoes instead of hand-copied numbers, PyCall.jl can call Echoes directly from a Julia session. This requires Echoes installed in some Python environment — any interpreter where import echoes succeeds will do, conda or otherwise:

julia

# One-time setup: point PyCall at the interpreter where `import echoes`
# succeeds, then rebuild the binding.
ENV["PYTHON"] = raw"/path/to/that/python"
using Pkg
Pkg.build("PyCall")

julia

using PyCall

py"""
from echoes import *

def run_echoes(phi, scheme_name):
    _SCHEMES = {"MT": MT, "SC": SC, "DIFF": DIFF}
    scheme = _SCHEMES[scheme_name]
    myrve = rve(matrix="SOLID")
    myrve["SOLID"] = ellipsoid(shape=spheroidal(1.), symmetrize=[ISO],
                                prop={"C": stiff_kmu(72., 32.)})
    myrve["PORE"] = ellipsoid(shape=spheroidal(1.), symmetrize=[ISO],
                               prop={"C": stiff_kmu(1e-6, 1e-6)})
    myrve["PORE"].fraction = phi
    myrve["SOLID"].fraction = 1. - phi
    C = homogenize(prop="C", rve=myrve, scheme=scheme,
                   epsrel=1e-6, maxnb=300, select_best=True)
    return max(C.k, 0.), max(C.mu, 0.)
"""

k_e, μ_e = py"run_echoes"(0.3, "MT")

Running this snippet against Echoes 1.0 returns (33.46058186790858, 17.626741619570563) at φ = 0.3, scheme = "MT", matching the reference table below to solver tolerance. This is the same pattern used throughout scripts/bench_echoes/ (see Production cross-checks) — a py"""...""" block defining a small helper, called from Julia like an ordinary function.

The benchmark result

The table below reproduces the captured Echoes 1.0 output for the porous benchmark above (ks, μs = 72, 32; φ = 0.0, 0.1, …, 1.0), obtained exactly as shown in the previous section. It is embedded here as literal data — reading it back requires no Echoes installation — so this page renders identically whether or not Echoes is available:

julia

# Captured from a run of Echoes 1.0 (see the run_echoes snippet above),
# once per (φ, scheme) pair.
echoes_ref = Dict(
    "MT" => (
        k = [72.0, 55.443851, 43.065421, 33.460582, 25.791046, 19.525425, 14.31056, 9.90258, 6.127661, 2.858562, 1.0e-6],
        μ = [32.0, 26.415585, 21.685158, 17.626742, 14.106633, 11.024391, 8.303101, 5.882864, 3.716342, 1.765626, 1.0e-6],
    ),
    "DIFF" => (
        k = [72.0, 54.60602, 40.616795, 29.414749, 20.528772, 13.593911, 8.323876, 4.491934, 1.917566, 0.457355, 1.0e-6],
        μ = [32.0, 26.168105, 20.864777, 16.11255, 11.933231, 8.347851, 5.376596, 3.038716, 1.352384, 0.334464, 1.0e-6],
    ),
    "SC" => (
        k = [72.0, 53.608383, 37.156454, 22.689248, 10.27169, 0.076007, 6.0e-6, 3.0e-6, 2.0e-6, 1.0e-6, 1.0e-6],
        μ = [32.0, 25.866249, 19.629147, 13.264365, 6.737918, 0.056937, 5.0e-6, 2.0e-6, 2.0e-6, 1.0e-6, 1.0e-6],
    ),
)

MFH, computed live, against those captured references:

julia

mfh_schemes = Dict(
    "MT" => MoriTanaka(),
    "DIFF" => DifferentialScheme(; nsteps = 300),
    "SC" => SelfConsistent(; reltol = 1.0e-6, abstol = 0.0, maxiters = 300, select_best = true),
)

_relerr(a, b) = abs(b) < 1.0e-9 ? abs(a - b) : abs(a - b) / abs(b)

plt = plot(;
    xlabel = "φ (porosity)", ylabel = "k_hom",
    legend = :topright, framestyle = :box, size = (760, 480),
)
for (name, color) in (("MT", :black), ("DIFF", :orange), ("SC", :red))
    k_jl = [max(k_mu(best_fit_iso(homogenize(build_rve(φ), mfh_schemes[name], :C)))[1], 0.0) for φ in φs]
    plot!(plt, φs, k_jl; label = "MFH $name", color = color, lw = 2)
    scatter!(plt, φs, echoes_ref[name].k; label = "Echoes $name", color = color, marker = :circle, markersize = 4)
end
plt

julia

for name in ("MT", "DIFF", "SC")
    scheme = mfh_schemes[name]
    println("== $name ==")
    for (i, φ) in enumerate(φs)
        k_jl, _ = k_mu(best_fit_iso(homogenize(build_rve(φ), scheme, :C)))
        k_jl = max(k_jl, 0.0)
        k_e = echoes_ref[name].k[i]
        println("  φ=", round(φ, digits = 2), "  k_MFH=", round(k_jl, digits = 6),
                "  k_Echoes=", k_e, "  relerr=", round(_relerr(k_jl, k_e), sigdigits = 3))
    end
end

== MT ==
  φ=0.0  k_MFH=72.0  k_Echoes=72.0  relerr=1.97e-16
  φ=0.1  k_MFH=55.443851  k_Echoes=55.443851  relerr=3.68e-9
  φ=0.2  k_MFH=43.065421  k_Echoes=43.065421  relerr=8.55e-9
  φ=0.3  k_MFH=33.460582  k_Echoes=33.460582  relerr=3.95e-9
  φ=0.4  k_MFH=25.791046  k_Echoes=25.791046  relerr=7.53e-9
  φ=0.5  k_MFH=19.525425  k_Echoes=19.525425  relerr=1.07e-8
  φ=0.6  k_MFH=14.31056  k_Echoes=14.31056  relerr=5.32e-9
  φ=0.7  k_MFH=9.90258  k_Echoes=9.90258  relerr=1.42e-8
  φ=0.8  k_MFH=6.127661  k_Echoes=6.127661  relerr=6.19e-8
  φ=0.9  k_MFH=2.858562  k_Echoes=2.858562  relerr=6.34e-8
  φ=1.0  k_MFH=1.0e-6  k_Echoes=1.0e-6  relerr=6.95e-9
== DIFF ==
  φ=0.0  k_MFH=72.0  k_Echoes=72.0  relerr=1.97e-16
  φ=0.1  k_MFH=54.611188  k_Echoes=54.60602  relerr=9.46e-5
  φ=0.2  k_MFH=40.631978  k_Echoes=40.616795  relerr=0.000374
  φ=0.3  k_MFH=29.43971  k_Echoes=29.414749  relerr=0.000849
  φ=0.4  k_MFH=20.56085  k_Echoes=20.528772  relerr=0.00156
  φ=0.5  k_MFH=13.62948  k_Echoes=13.593911  relerr=0.00262
  φ=0.6  k_MFH=8.359112  k_Echoes=8.323876  relerr=0.00423
  φ=0.7  k_MFH=4.523194  k_Echoes=4.491934  relerr=0.00696
  φ=0.8  k_MFH=1.94151  k_Echoes=1.917566  relerr=0.0125
  φ=0.9  k_MFH=0.470894  k_Echoes=0.457355  relerr=0.0296
  φ=1.0  k_MFH=1.0e-6  k_Echoes=1.0e-6  relerr=1.66e-10
== SC ==
  φ=0.0  k_MFH=72.0  k_Echoes=72.0  relerr=1.97e-16
  φ=0.1  k_MFH=53.608383  k_Echoes=53.608383  relerr=3.58e-10
  φ=0.2  k_MFH=37.156446  k_Echoes=37.156454  relerr=2.13e-7
  φ=0.3  k_MFH=22.689245  k_Echoes=22.689248  relerr=1.47e-7
  φ=0.4  k_MFH=10.271684  k_Echoes=10.27169  relerr=6.05e-7
  φ=0.5  k_MFH=0.076512  k_Echoes=0.076007  relerr=0.00665
  φ=0.6  k_MFH=6.0e-6  k_Echoes=6.0e-6  relerr=0.00386
  φ=0.7  k_MFH=3.0e-6  k_Echoes=3.0e-6  relerr=0.0575
  φ=0.8  k_MFH=2.0e-6  k_Echoes=2.0e-6  relerr=0.101
  φ=0.9  k_MFH=1.0e-6  k_Echoes=1.0e-6  relerr=0.35
  φ=1.0  k_MFH=1.0e-6  k_Echoes=1.0e-6  relerr=0.0

Scheme	Agreement	Source of the residual
`Mori-Tanaka`		floating-point floor; closed form on both sides
`Differential`	sub-percent to , a few percent beyond	independent step counts and quadrature on the same ODE
`SelfConsistent`	close to	past percolation both sides sit at solver tolerance (), so the relative error stops being meaningful

Production cross-checks

scripts/bench_echoes/ runs the same PyCall pattern systematically; tolerances and measured agreement are in Cross-validation, timings in Performance vs Echoes.

Two tutorials carry a captured Echoes reference in the page itself, in the literal-data style used above, and are worth reading next because each one turns up a place where the two libraries do not line up trivially:

Page	What the cross-check shows
Crack distributions: isotropic or parallel	`MoriTanaka` and `SelfConsistent` reproduce Echoes to on both an isotropic and a single-orientation crack population — but `AsymmetricSelfConsistent` is a different fixed point, not the translation of Echoes' `SC`, and it percolates at about half the crack density
Ageing creep: loading age against inclusion shape	the ageing-viscoelastic counterpart of the table above: `ViscoLaw(J, :creep)` against Echoes' `visco_prop={"C": (Js, CREEP)}`, swept over loading ages and aspect ratios

Without Echoes or PyCall, each cross-check has a PyCall-free counterpart: committed *_python.json dumps, and pure-Julia transcriptions such as scripts/55_ageing_creep_dirichlet_chains.jl.

From Echoes to MeanFieldHomogenization ​

API translation table ​

Same problem, both sides ​

Calling Echoes from Julia with PyCall ​

The benchmark result ​

Production cross-checks ​

From Echoes to MeanFieldHomogenization

API translation table

Same problem, both sides

Calling Echoes from Julia with PyCall

The benchmark result

Production cross-checks