API — Laplace-Carson viscoelasticity

The non-ageing half of MeanFieldHomogenization.Viscoelasticity: numerical Laplace inversion, the rheological model catalog, the exact Kelvin ⇄ Maxwell conversion, and the homogenization driver that ties them together.

See the theory, the model manual and the inversion manual.

Numerical inversion

MeanFieldHomogenization.Viscoelasticity.AbstractLaplaceInversion Type

AbstractLaplaceInversion

Root supertype for numerical inverse-Laplace algorithms: GaverStehfest, FixedTalbot, TalbotTrefethen and DeHoog.

All are consumed through inverse_laplace and inverse_carson.

MeanFieldHomogenization.Viscoelasticity.GaverStehfest Type

GaverStehfest(N::Int = 16)

Gaver-Stehfest inversion,

with the Salzer weights .

All nodes are real and positive. That is why this is the default: a transform defined only for real p works, an entire homogenization scheme run through homogenize_lc stays in real arithmetic, and ForwardDiff sees plain Dual numbers rather than Complex{Dual}.

N must be even with 4 ≤ N ≤ 40.

More terms is not more accurate

The weights alternate in sign with magnitudes up to times the answer, so the method needs roughly 2.3 M digits of working precision to return M correct ones. In Float64 that puts a hard ceiling on N:

N	relative error on   at t = 1, in Float64
8	8e-3
12	4e-4
16 (default)	1e-5
18	2e-6 — the optimum
22	2e-3
26	4e-1
30	1e+3 — and it keeps growing

Past N ≈ 18 the round-off grows — the single most surprising property of the method, and the reason N is validated rather than left free.

Working in BigFloat removes the ceiling entirely: the weights are cached as exact Rational{BigInt} and converted to the working type on demand, so under setprecision(256) the error keeps falling with N all the way to 1e-14 at N = 34 — where Float64 is at 1e5. See the inversion tutorial §2 for the plot.

At its best on completely monotone kernels — creep and relaxation functions — and poor on oscillatory ones: on at N = 16 the relative error is O(1), where FixedTalbot reaches 1e-12.

Measured on    with τ = 1, N = 16:

t/τ	relative error, E_∞ > 0	relative error, E_∞ = 0
0.01	4e-8	6e-8
1	3e-9	5e-7
10	3e-5	0.4
40	3e-6	5e11

The second column is not a defect of this method — see the note on the tail under inverse_laplace. The first is the regime viscoelastic solids actually live in.

MeanFieldHomogenization.Viscoelasticity.FixedTalbot Type

FixedTalbot(N::Int = 24)

Fixed-Talbot inversion (Abate & Valkó): the Bromwich contour is deformed into the cotangent curve

on which decays fast enough for the midpoint trapezoidal rule to converge geometrically. Near machine precision with N = 24 evaluations on meromorphic transforms — Prony series, Zener, Burgers.

This is DEFAULT_INVERSION, and on every kernel tested it is also the most accurate of the four — around 1e-12–1e-13 relative, for 24 evaluations of the transform:

transform	FixedTalbot(24)	GaverStehfest(16)
1/(p+a) (exponential)	3e-12	1e-5
1/p³ (polynomial)	1e-14	3e-7
1/√p (branch cut)	4e-12	4e-7
1/(p²+ω²) (oscillatory)	1e-12	O(1)
2S2P1D creep kernel	5e-13	3e-8

Branch cuts are fine — that is what the contour is for

The cotangent curve is a Hankel-type contour: it wraps around the negative real axis rather than crossing it, which is precisely how Talbot quadratures are designed to handle singularities there. Fractional models (ScottBlair, HuetSayegh, Model2S2P1D, Rabotnov) invert at full accuracy — verified against their exact   pairs in test/Viscoelasticity/.

What it cannot do is reach a singularity in the right half-plane; use TalbotTrefethen(; shift = …) for that.

Requires the transform to accept complex arguments. When it must stay real — so that a whole homogenization scheme runs in real arithmetic — use GaverStehfest instead.

MeanFieldHomogenization.Viscoelasticity.TalbotTrefethen Type

TalbotTrefethen(N::Int = 24; shift = 0.0)

Talbot inversion on the Trefethen-Weideman-Schmelzer contour,

with  ,  ,  ,  . Converges as , and like FixedTalbot it handles branch cuts on   at full accuracy. This is the variant used by the ECHOES Python reference (tests/python/creep/laplace_inversion.py), so results can be cross-checked against it directly.

It differs from FixedTalbot in one respect that matters: its contour is tuned for decaying kernels and degrades badly on oscillatory ones (6e-2 on at t = 2, against 9e-13 for FixedTalbot). Prefer FixedTalbot unless cross-checking, or unless shift is needed.

shift moves the contour to the right when the transform has poles on the positive real axis.

MeanFieldHomogenization.Viscoelasticity.DeHoog Type

DeHoog(; N = 16, T = nothing, tol = 1e-9)

De Hoog-Knight-Stokes inversion: the Bromwich integral is discretized into a Fourier series on the vertical line  ,

whose slow convergence is accelerated by the quotient-difference algorithm and summed as a continued fraction. Ten to thirteen digits in practice.

Because the line stays in the right half-plane it never crosses a branch cut on   — unlike the Talbot family — so this is the general-purpose choice for fractional models when more accuracy than GaverStehfest is wanted.

T is the scaling period. The nodes depend on T but not on t, so one node set — one pass over F — can serve several times at once. What limits the sharing is that the relative accuracy depends on the ratio t/T alone:

t/T	relative error (N = 16, tol = 1e-9)
0.5 (i.e. T = 2t)	≈ 1e-9
0.15	≈ 2e-10
0.05	≈ 1e-6
0.005	≈ 1e-2
0.0005	≈ 0.4 — meaningless

and raising N barely helps below t/T ≈ 0.05. A single node set therefore covers a window of times spanning a factor of about three, not several decades.

Two modes follow:

• T = nothing (the default) — on a grid, the times are sorted and split into blocks spanning at most a factor of three, each block getting its own T = 2 t_max and its own single pass over F. Accuracy is uniform, and a 200-point grid over seven decades costs roughly 15 × (2N+1) ≈ 500 evaluations of F instead of 200 × (2N+1) ≈ 6600. That is the reason to reach for DeHoog when each evaluation of F is a homogenization scheme.

• T given explicitly — one node set for everything, 2N + 1 evaluations total. Only do this when the grid really is narrow; a warning is emitted for any time falling below t/T = 0.15.

Not a weighted sum

Unlike the other three, de Hoog is not linear in the transform values: the quotient-difference tables divide consecutive samples by one another. Tensor- and matrix-valued transforms are therefore inverted component by component, and a component that vanishes identically is short-circuited to zero rather than run through a 0/0.

MeanFieldHomogenization.Viscoelasticity.DEFAULT_INVERSION Constant

DEFAULT_INVERSION

The inversion method used when none is given: FixedTalbot(24).

Chosen on measured accuracy: it is the only one of the four that stays near 1e-12 on every kernel tested — exponential, polynomial, branch-cut, oscillatory and the fractional 2S2P1D pair alike (see FixedTalbot for the table).

Reasons to override it:

• GaverStehfest when the transform must be evaluated at real p only — that keeps a whole homogenize_lc sweep in real arithmetic, costs a third fewer evaluations per point, and is the one way to use SelfConsistent(algorithm = NewtonDefault()) in the Laplace-Carson route;

• DeHoog when inverting on a grid whose every point costs a homogenization, since one node set serves a whole block of times;

• TalbotTrefethen to cross-check against the ECHOES reference, or when the transform has poles in the right half-plane.

MeanFieldHomogenization.Viscoelasticity.inverse_laplace Function

inverse_laplace(F, t, method = DEFAULT_INVERSION)
inverse_laplace(F, times::AbstractVector, method = DEFAULT_INVERSION)

Numerically invert the Laplace transform F, returning f(t) with

F is any callable p -> F(p) whose value supports + and multiplication by a scalar. Scalars, TensND tensors of any symmetry class and 6×6 Mandel matrices all work, and the symmetry class of the result is the class F returns — a TensTI transform gives back a TensTI with its axis preserved.

method selects the algorithm; see AbstractLaplaceInversion. FixedTalbot, TalbotTrefethen and DeHoog call F with Complex arguments; GaverStehfest only with real positive ones.

The vector form returns a Vector of results. With DeHoog(; T = ...) it shares one node set across the whole grid, which matters when each evaluation of F is a homogenization scheme.

Throws DomainError for t ≤ 0: every method places its nodes at O(1/t).

Name clash with Symbolics.jl

Symbolics exports a function also called inverse_laplace — a five-argument symbolic transform, inverse_laplace(expr, F, s, f, t), used for solving ODEs. The two are different functions, so using MeanFieldHomogenization, Symbolics makes the bare name ambiguous and Julia refuses it. Qualify whichever you mean:

MeanFieldHomogenization.inverse_laplace(F, t)     # this one
Symbolics.inverse_laplace(expr, F, s, f, t)       # theirs

inverse_carson — the one to use for anything in the rheology catalog — is not affected.

Accuracy is absolute, not relative, in the tail

Every algorithm here controls the error against the scale of f, not against f(t) at the point asked for. Once f has decayed many orders of magnitude below its initial value the relative error is unbounded: on at t = 40, where the function is 1e-17, FixedTalbot(24) is off by a relative 1e4 — an absolute 1e-13, which is exactly what it promises.

This is intrinsic to numerical Laplace inversion, not a defect of one method, and it is why viscoelastic solids — whose relaxation function settles on a plateau   rather than decaying to zero — are the comfortable case. For a fluid, read the tail as an absolute quantity.

Examples

julia> using MeanFieldHomogenization

julia> isapprox(inverse_laplace(p -> 1 / (p + 2), 1.0), exp(-2.0); rtol = 1e-6)
true

See also inverse_carson, inverse_carson_rate.

MeanFieldHomogenization.Viscoelasticity.inverse_carson Function

inverse_carson(Fstar, t, method = DEFAULT_INVERSION)
inverse_carson(Fstar, times::AbstractVector, method = DEFAULT_INVERSION)

Invert the Laplace-Carson transform

i.e. inverse_laplace(p -> Fstar(p) / p, t, method).

Laplace-Carson is the transform of choice in viscoelasticity because it maps a constant to itself: an elastic modulus is its own transform, so the correspondence principle reads C*(p) in place of C with no extra factor. Every model in the rheology catalog exposes its transform in this convention through carson_relaxation and carson_creep.

Examples

julia> using MeanFieldHomogenization

julia> R = zener_maxwell(2.0, 3.0, 1.0);      # E_∞ = 2, E₁ = 3, τ₁ = 1

julia> r = inverse_carson(p -> carson_relaxation(R, p), 0.5);

julia> isapprox(r, relaxation(R, 0.5); rtol = 1e-6)
true

MeanFieldHomogenization.Viscoelasticity.inverse_carson_rate Function

inverse_carson_rate(Fstar, t, method = DEFAULT_INVERSION; f_glassy)

Return the time derivative ḟ(t) of the function whose Laplace-Carson transform is Fstar, without differentiating the inversion.

Because L{ḟ}(p) = p f̂(p) - f(0⁺) = f*(p) - f(0⁺), the rate is itself an ordinary inverse Laplace transform:

f_glassy is f(0⁺) = lim_{p→∞} f*(p), available in closed form for every model in the catalog via glassy_modulus.

Relation to ForwardDiff.derivative(t -> ..., t)

Differentiating the inversion with respect to t also works — every method here accepts a Dual time — and on the kernels tested the two routes are of comparable accuracy. The identity is worth having anyway because it asks nothing of the transform beyond what it already provides: no Dual is pushed into the nodes, so it applies when F is not differentiable, and it is the cheaper route when the rate is wanted alongside the value.

Internals

MeanFieldHomogenization.Viscoelasticity._accumulate Function

_accumulate(values, weights)

Form Σ_k weights[k] * values[k], seeding the accumulator with the first term rather than with zero(...).

That single choice carries two unrelated correctness requirements at once:

• ForwardDiff — a zero(T) seed derived from the node type would pin the accumulator to a plain float and silently drop the partials of the transform values;

• TensND — Base.zero(::AbstractTens{4,dim,T}) returns a TensISO whatever the input class is (the @eval loop over one/zero in TensND/src/tens_isotropic.jl). Seeding with zero(A) would therefore collapse a TensTI transform onto the isotropic class and lose its axis. Never calling zero sidesteps the trap entirely.

MeanFieldHomogenization.Viscoelasticity._realpart Function

_realpart(x)

Discard the imaginary part left over by a numerical inversion, whatever the container. Every inversion method here evaluates the transform on complex nodes, so a real-valued original comes back with an imaginary residue of the order of the quadrature error; the caller wants the real object back, not a complex one.

The tensor methods rebuild through TensND._rebuild rather than through get_array, which is what keeps a TensISO/TensTI/TensOrtho in its own class with its axis or frame intact — dropping to a generic Tens would lose the symmetry the scheme just established.

MeanFieldHomogenization.Viscoelasticity._decompose Function

_decompose(x) -> (components::Vector, rebuild)

Split a transform value into a vector of scalar components together with a closure putting them back into the same shape and symmetry class.

Used only by DeHoog, whose quotient-difference acceleration divides samples by one another and is therefore not linear in the transform values. The other three methods never need this: they are weighted sums, so they run on the value type directly.

MeanFieldHomogenization.Viscoelasticity._gs_exact_weights Function

_gs_exact_weights(N) -> Vector{Rational{BigInt}}

Salzer weights of the Gaver-Stehfest scheme, exactly,

MeanFieldHomogenization.Viscoelasticity._scalar_float Function

_scalar_float(T)

The plain floating-point type underlying T, peeling ForwardDiff.Dual wrappers. The quadrature weights are built in this type: they never carry partials (they depend on N and t only), and promotion does the rest.

MeanFieldHomogenization.Viscoelasticity._plain_value Function

_plain_value(x)

Strip every ForwardDiff.Dual layer, returning the underlying float.

Used only to choose DeHoog's scaling period T. T is a free algorithmic parameter — it selects the node set, not the mathematics — so freezing it at the value level is legitimate and keeps the Dual confined to the one place where t enters analytically, exp(iπt/T). Letting T carry partials would instead push them into S(Tscale) and fail.

MeanFieldHomogenization.Viscoelasticity._dehoog_from_values Function

_dehoog_from_values(values, t, method, Tscale)

De Hoog-Knight-Stokes evaluation from transform samples already taken on the 2N+1 nodes.

The Fourier series — with  ,   and   — is turned into a continued fraction by the quotient-difference algorithm, summed by the standard three-term recurrence, and closed with the remainder estimate R2M that removes the leading truncation error.

The QD tables divide samples by one another, so this is not linear in the transform values: tensor- and matrix-valued transforms are inverted component by component through _decompose.

MeanFieldHomogenization.Viscoelasticity._dehoog_blocks Function

_dehoog_blocks(order, times, ratio) -> Vector{UnitRange}

Split the sorted indices order into consecutive runs whose times span at most ratio, so each run can share one de Hoog node set.

MeanFieldHomogenization.Viscoelasticity._DEHOOG_MIN_RATIO Constant

_DEHOOG_MIN_RATIO

Smallest t/T at which DeHoog still returns close to full accuracy. Measured on exponential and power-law kernels: the relative error is ≈ 2e-10 at t/T = 0.15 and ≈ 1e-6 already at t/T = 0.05, essentially independent of N.

A shared node set therefore spans a factor 0.5 / 0.15 ≈ 3.3 in time, which is what _dehoog_blocks uses.

The model interface

MeanFieldHomogenization.Viscoelasticity.AbstractRheology Type

AbstractRheology

Root supertype for scalar linear viscoelastic models — a spring, a Prony series, 2S2P1D, and everything in between.

A concrete model must define carson_relaxation; it should also define glassy_modulus and equilibrium_modulus, which are limits of that transform and are almost always elementary. Anything else it knows in closed form (relaxation, creep, carson_creep) is an override; what it does not know is obtained by numerical inversion.

Lift a pair of scalar models to a fourth-order isotropic tensor with iso_rheology / iso_rheology_E_nu; the result is an AbstractTensorRheology, which answers the same generics with tensors and additionally feeds the ageing time-domain pipeline through ViscoLaw(model).

See the rheological model catalog.

MeanFieldHomogenization.Viscoelasticity.AbstractTensorRheology Type

AbstractTensorRheology

Root supertype for tensor-valued viscoelastic models — a scalar model in each independent symmetry channel. The only concrete subtype at present is IsoRheology.

MeanFieldHomogenization.Viscoelasticity.carson_relaxation Function

carson_relaxation(model, p)

The Laplace-Carson relaxation transform    — the one method a model must implement.

p may be real (as GaverStehfest needs), complex (p = iω gives the complex modulus), or a ForwardDiff.Dual, so implementations must stay generic in the argument type.

For an AbstractTensorRheology the value is a fourth-order tensor.

MeanFieldHomogenization.Viscoelasticity.carson_creep Function

carson_creep(model, p)

The Laplace-Carson creep transform .

In the Carson domain creep and relaxation are exact reciprocals,   , which is the default implementation. Models whose creep transform is the simpler of the two (PronyCreep, burgers) override it and let carson_relaxation be the derived one.

MeanFieldHomogenization.Viscoelasticity.relaxation Function

relaxation(model, t)

The relaxation function : the stress response to a unit strain step applied at t = 0.

Falls back to numerically inverting carson_relaxation with default_inversion; models with a closed form override it. t = 0 is answered analytically by glassy_modulus — the inversion itself has a pole there.

MeanFieldHomogenization.Viscoelasticity.creep Function

creep(model, t)

The creep compliance : the strain response to a unit stress step applied at t = 0.

Falls back to numerically inverting carson_creep; models with a closed form override it.

MeanFieldHomogenization.Viscoelasticity.complex_modulus Function

complex_modulus(model, ω)

The complex modulus  , the quantity a dynamic-mechanical test measures.

Its modulus abs(E*) is the norm of the complex modulus and its argument angle(E*) the phase angle; storage_modulus, loss_modulus and loss_factor name the usual derived quantities.

MeanFieldHomogenization.Viscoelasticity.storage_modulus Function

storage_modulus(model, ω)

  — the in-phase, elastically stored part.

MeanFieldHomogenization.Viscoelasticity.loss_modulus Function

loss_modulus(model, ω)

  — the out-of-phase, dissipated part.

MeanFieldHomogenization.Viscoelasticity.loss_factor Function

loss_factor(model, ω)

  — the phase angle's tangent. Non-negative for any passive material, which makes it a cheap sanity check on a fitted spectrum.

MeanFieldHomogenization.Viscoelasticity.glassy_modulus Function

glassy_modulus(model)

  — the instantaneous ("glassy") modulus.

Closed form for every model in the catalog. Besides being the t = 0 value of relaxation, it is the f_glassy argument of inverse_carson_rate.

MeanFieldHomogenization.Viscoelasticity.equilibrium_modulus Function

equilibrium_modulus(model)

  — the relaxed ("equilibrium" or "static") modulus.

Zero exactly when the model is a fluid: stress relaxes away completely and the creep compliance grows without bound. See is_fluid.

MeanFieldHomogenization.Viscoelasticity.is_fluid Function

is_fluid(model) -> Bool

true when equilibrium_modulus vanishes, i.e. the model contains a dashpot in series with everything else and creeps indefinitely.

The distinction is structural rather than cosmetic: a fluid's creep compliance carries a term linear in t, which in the Carson domain is the pole φ/p of PronyCreep, and it is what decides whether kelvin_to_maxwell and maxwell_to_kelvin look for a root in the outermost interval.

MeanFieldHomogenization.Viscoelasticity.default_inversion Function

default_inversion(model) -> AbstractLaplaceInversion

The inversion method relaxation and creep use when the model has no closed form and the caller gave none.

Defaults to DEFAULT_INVERSION; a model overrides it when its transform has a feature one algorithm handles better.

Discrete spectra, and the exact conversion

MeanFieldHomogenization.Viscoelasticity.PronyRelaxation Type

PronyRelaxation(E_inf, E, tau)

A generalized Maxwell chain: an equilibrium spring E_inf in parallel with branches (E[i], tau[i]), each a spring and a dashpot in series.

E_inf == 0 means a fluid: the stress relaxes away completely.

Why the degenerate branches are separate fields

The ECHOES Python reference stores the equilibrium term inside the spectrum as a branch with τ = Inf (and, in the dual type, τ = 0). That forces Inf into a vector which is then sorted, exponentiated and differentiated — and Inf is poison for ForwardDiff, where Inf * 0 gives NaN in the partials. Keeping E_inf as its own scalar removes every special case.

Construction sorts the branches by increasing tau and merges any that agree to within merge_tol (relative), summing their moduli — two branches with the same relaxation time are one branch, not an error.

See maxwell_to_kelvin for the conversion, prony_fit_relaxation to obtain one by fitting an arbitrary transform, and zener_maxwell for the one-branch case.

MeanFieldHomogenization.Viscoelasticity.PronyCreep Type

PronyCreep(J_0, J, tau, phi = zero(J_0))

A generalized Kelvin chain: an instantaneous spring of compliance J_0 in series with branches (J[i], tau[i]), each a spring and a dashpot in parallel, and — when phi > 0 — a dashpot of fluidity phi.

Fluidity, not viscosity

phi = 1/η rather than η, so a solid is phi = 0 rather than η = Inf. Zero is differentiable and sorts and promotes like any other number; Inf does none of those things. It also makes the two Prony types exactly symmetric:

phi == 0    ⟺  solid  ⟺  E_inf > 0
E_inf == 0  ⟺  fluid  ⟺  phi > 0

See kelvin_to_maxwell, prony_fit_creep, burgers and zener_kelvin.

MeanFieldHomogenization.Viscoelasticity.PRONY_MERGE_TOL Constant

PRONY_MERGE_TOL

Relative spacing below which two branches of a Prony spectrum are considered the same branch and merged on construction. Also the scale below which maxwell_to_kelvin warns that the converted spectrum, though still mathematically exact, is ill-conditioned.

MeanFieldHomogenization.Viscoelasticity.PRONY_FLUID_TOL Constant

PRONY_FLUID_TOL

Relative size below which an equilibrium modulus (or a fluidity) counts as zero, so the model is treated as a fluid (resp. a solid). Compared against the glassy modulus, so it is scale-free.

MeanFieldHomogenization.Viscoelasticity.maxwell_to_kelvin Function

maxwell_to_kelvin(m::PronyRelaxation; rtol = 1e-14) -> PronyCreep

Convert a generalized Maxwell chain into the exactly equivalent generalized Kelvin chain, so that J*(p) R*(p) = 1 identically.

How

In σ = -1/p, the transform becomes    , which is strictly increasing between consecutive poles and runs from -∞ to +∞ across every gap. The retardation times are its zeros, and they interlace the relaxation times:

     0 <  τ₁  <  σ₁  <  τ₂  <  σ₂  < … <  τ_m  < [σ_m]

with the last root present exactly when E_∞ > 0. There is deliberately no root in (0, τ₁): Φ(0⁺) = E_glassy > 0 and Φ only increases from there — the bracketing is not symmetric with kelvin_to_maxwell, which does have a root below its first pole.

Each root is thus isolated before any arithmetic happens, and bisection in log σ finds it whatever the number of branches or their spread. The residues follow from R*'(-1/σ) = σ²Φ'(σ):

and Φ' > 0 makes every J_j positive by construction rather than by luck.

Fluids

When E_∞ = 0 the outermost root recedes to infinity: that is the series dashpot. It is picked up exactly, as   , and the result has one branch fewer than the input.

Differentiability

ForwardDiff traverses the conversion: the roots are lifted by the implicit function theorem (see _root_with_ad), so gradients of J_j and τ^K_j with respect to (E_∞, E, τ) are available.

Examples

julia> using MeanFieldHomogenization

julia> m = PronyRelaxation(8.0, [3.0, 17.0], [12.0, 23.0]);

julia> k = maxwell_to_kelvin(m);

julia> isapprox(carson_creep(k, 0.37) * carson_relaxation(m, 0.37), 1.0; rtol = 1e-12)
true

See also kelvin_to_maxwell, prony_fit_creep.

MeanFieldHomogenization.Viscoelasticity.kelvin_to_maxwell Function

kelvin_to_maxwell(k::PronyCreep; rtol = 1e-14) -> PronyRelaxation

Convert a generalized Kelvin chain into the exactly equivalent generalized Maxwell chain, so that R*(p) J*(p) = 1 identically. The inverse of maxwell_to_kelvin, and the operation Kelvin2Maxwell.py performs symbolically in ECHOES.

How

In σ = -1/p,      is strictly decreasing between consecutive poles. Its zeros are the relaxation times, and here — unlike the other direction — there is one below the first pole, because Ψ(0⁺) = J_0 > 0 while Ψ(τ₁⁻) = -∞:

     0 <  σ₁  <  τ₁  <  σ₂  <  τ₂  < … <  τ_n  < [σ_{n+1}]

the last root existing exactly when φ > 0. So a solid Kelvin chain of n branches gives n Maxwell branches with E_∞ > 0, and a fluid one gives n+1 branches with E_∞ = 0 — which is the degree count of the rational transform, as it should be.

The residues, from J*'(-1/σ) = σ²Ψ'(σ):

with Ψ' < 0 making every E_j positive by construction.

For a solid the result satisfies a free consistency identity — both sides are R*(0):

which is checked internally and would catch any sign slip in the residues.

Examples

The Burgers model is a fluid Kelvin chain with one branch, so it must convert to a two-branch Maxwell chain with no equilibrium spring:

julia> using MeanFieldHomogenization

julia> b = burgers(1.0, 3.0, 2.0, 6.0);

julia> r = kelvin_to_maxwell(b);

julia> length(r), isapprox(equilibrium_modulus(r), 0.0; atol = 1e-12)
(2, true)

See also maxwell_to_kelvin, prony_fit_relaxation.

MeanFieldHomogenization.Viscoelasticity.prony_fit_relaxation Function

prony_fit_relaxation(Rstar, taus; points = nothing, E_inf = nothing,
                     nonneg = true, lambda = 0.0) -> PronyRelaxation

Fit a generalized Maxwell chain with the prescribed relaxation times taus to an arbitrary Laplace-Carson relaxation transform Rstar, by least squares on a set of collocation points.

The model is linear in its unknowns,

so the fit is one linear solve — no iteration, no starting guess.

• points — the Carson variables to collocate at. Defaults to 2 * length(taus) values log-spaced over [1/(10 τ_max), 10/τ_min], which brackets the spectrum by a decade on each side.

• E_inf — the equilibrium modulus. Defaults to Rstar evaluated at 1/(100 τ_max).

  Note

  The ECHOES reference uses rstar(1e-100) for this. That overflows p^{-k} for any fractional model — 2S2P1D included — so the default here is tied to the spectrum's own slowest time instead.

• nonneg — constrain X_i ≥ 0 (the default). A non-negative spectrum is what makes the fitted function completely monotone, hence passive; an unconstrained fit routinely produces small negative moduli that give a negative loss_factor somewhere.

• lambda — Tikhonov weight, useful when the times are packed more densely than the data can resolve.

Not differentiable

With nonneg = true the solve is an active-set method: the returned coefficients are a piecewise-smooth function of the input with a combinatorial switch in the middle. Treat a fit as a calibration step that produces a model, and differentiate the model — not the fit. An unconstrained fit (nonneg = false) is a plain linear solve and is differentiable.

See prony_fit_creep for the dual, and maxwell_to_kelvin for the exact conversion once a chain is in hand.

MeanFieldHomogenization.Viscoelasticity.prony_fit_creep Function

prony_fit_creep(Jstar, taus; points = nothing, J_0 = nothing, phi = 0.0,
                nonneg = true, lambda = 0.0) -> PronyCreep

The creep-side twin of prony_fit_relaxation: fit a generalized Kelvin chain to a Laplace-Carson creep transform Jstar,

J_0 defaults to Jstar at 100/τ_min, i.e. well above the fastest retardation time.

phi is not fitted: a φ/p pole is a qualitatively different object from the rest of the sum, and reading it off noisy samples is unreliable. Pass it explicitly when the material is a fluid — it is 1/η of the series dashpot, and for a transform known in closed form it is lim_{p→0} p J*(p).

Internals

MeanFieldHomogenization.Viscoelasticity._phi_maxwell Function

_phi_maxwell(m::PronyRelaxation, σ)

     .

Its zeros are the retardation times of the equivalent Kelvin chain. Strictly increasing between consecutive poles, with Φ(0⁺) = E_glassy > 0 and Φ(∞) = E_∞.

MeanFieldHomogenization.Viscoelasticity._dphi_maxwell Function

_dphi_maxwell(m::PronyRelaxation, σ)

  , strictly positive for a positive spectrum — which is why the converted compliances J_j = 1/(σ_j Φ'(σ_j)) come out positive by construction.

MeanFieldHomogenization.Viscoelasticity._psi_kelvin Function

_psi_kelvin(m::PronyCreep, σ)

      .

Its zeros are the relaxation times of the equivalent Maxwell chain. Strictly decreasing between consecutive poles, with Ψ(0⁺) = J_0 > 0.

MeanFieldHomogenization.Viscoelasticity._dpsi_kelvin Function

_dpsi_kelvin(m::PronyCreep, σ)

    , strictly negative — hence E_j = -1/(σ_j Ψ'(σ_j)) > 0.

MeanFieldHomogenization.Viscoelasticity._bracketed_root Function

_bracketed_root(f, lo, hi; rtol = 1e-14, maxiter = 200) -> Float64

Bisection in log σ on a bracket known to contain exactly one simple root.

Logarithmic rather than linear because relaxation spectra span six to ten decades: the geometric midpoint √(lo·hi) halves the relative interval, so the iteration count depends on the requested relative accuracy and not at all on where in the spectrum the root sits.

The caller guarantees the sign change; a violation is a bug in the bracketing tables, not a user error, so it throws.

MeanFieldHomogenization.Viscoelasticity._expand_bracket Function

_expand_bracket(f, lo, growth, maxsteps) -> hi

Find an upper end for a root known to lie somewhere in (lo, ∞) by repeatedly multiplying by growth until the sign flips. Used for the outermost interval, which is unbounded whenever the model has a degenerate branch at the slow end.

MeanFieldHomogenization.Viscoelasticity._shrink_bracket Function

_shrink_bracket(f, hi, shrink, maxsteps) -> lo

The mirror image of _expand_bracket for a root in (0, hi).

MeanFieldHomogenization.Viscoelasticity._root_with_ad Function

_root_with_ad(f, df, lo, hi) -> σ

Locate the root of f in [lo, hi] and hand back a value that carries the correct ForwardDiff partials.

The bisection runs on values only — an iteration count is discrete, and differentiating through it is meaningless. A single Newton step taken in the full Dual type then supplies the derivatives: since f(σ_val) ≈ 0 at the value level, σ - f(σ)/f'(σ) leaves the value alone and sets the partials to

which is exactly the implicit function theorem for a simple root. No nested duals, and no differentiation of the solver.

MeanFieldHomogenization.Viscoelasticity._sort_and_merge Function

_sort_and_merge(w, tau, merge_tol, who) -> (w, tau)

Sort a spectrum by increasing time and fold together branches whose times agree to within merge_tol relative, summing their weights.

Merging is the correct preprocessing rather than an error: two branches with the same characteristic time are indistinguishable, and their sum is the same material. It also guarantees the strict interlacing the conversion relies on — maxwell_to_kelvin needs τ_i < τ_{i+1} with no ties, or the brackets it bisects in would be empty.

MeanFieldHomogenization.Viscoelasticity._default_collocation Function

_default_collocation(taus) -> Vector{Float64}

2 length(taus) Carson variables log-spaced over [1/(10 τ_max), 10/τ_min].

Real and positive, so the transform is only ever sampled where every model in the catalog is defined — including the fractional ones, whose p^{-k} needs a branch choice off the positive real axis.

MeanFieldHomogenization.Viscoelasticity._lsq Function

_lsq(A, b, nonneg, lambda) -> x

Least-squares solve of A x ≈ b, optionally Tikhonov-regularized by lambda and optionally constrained to x ≥ 0.

MeanFieldHomogenization.Viscoelasticity._nnls Function

_nnls(A, b; maxiter, tol) -> x

Lawson-Hanson non-negative least squares: minimize ‖Ax - b‖₂ subject to x ≥ 0.

Implemented here rather than pulled in as a dependency — it is forty lines, and the alternative in the ECHOES reference is an nlopt call. Being an active-set method it is not differentiable; see the warning on prony_fit_relaxation.

The model catalog

Elementary elements

MeanFieldHomogenization.Viscoelasticity.Spring Type

Spring(E)

A Hookean spring: R(t) = E, J(t) = 1/E, R*(p) = E.

The elastic limit of every other model, and the natural way to give a phase a constant modulus in one channel while another channel relaxes — the bulk modulus of a bituminous binder, for instance.

Being its own Laplace-Carson transform is not an accident: it is the property that makes the Carson convention the natural one in viscoelasticity.

MeanFieldHomogenization.Viscoelasticity.Dashpot Type

Dashpot(η)

A Newtonian dashpot: J(t) = t/η, R*(p) = pη.

Its relaxation function is η δ(t), a distribution rather than a function, so relaxation throws for this model — as does glassy_modulus, which would be infinite. A dashpot is meant to be assembled with springs (MaxwellUnit, burgers), not loaded on its own.

MeanFieldHomogenization.Viscoelasticity.MaxwellUnit Type

MaxwellUnit(E, η)

A spring and a dashpot in series:

A fluid: the stress relaxes to zero and the strain grows without bound. Equivalent to PronyRelaxation(0, [E], [η/E]), and that is what kelvin_to_maxwell and friends produce; this type exists so the physical parameters (E, η) can be given directly.

MeanFieldHomogenization.Viscoelasticity.KelvinUnit Type

KelvinUnit(E, η)

A spring and a dashpot in parallel:

The strain is bounded but the instantaneous response is rigid, so the relaxation function again carries a Dirac and relaxation throws. Put it in series with a spring — zener_kelvin — for a usable solid.

MeanFieldHomogenization.Viscoelasticity.relaxation_time Function

relaxation_time(m::MaxwellUnit) -> Number

The time constant τ = η / E of a Maxwell unit: the time over which its stress decays to 1/e of its initial value under a held strain. It is the sole shape parameter of the unit's relaxation function, R(t) = E exp(-t/τ).

See retardation_time for the Kelvin counterpart.

MeanFieldHomogenization.Viscoelasticity.retardation_time Function

retardation_time(m::KelvinUnit) -> Number

The time constant τ = η / E of a Kelvin unit: the time over which its strain reaches 1 - 1/e of the asymptotic value under a held stress. It is the sole shape parameter of the unit's creep function, J(t) = (1 - exp(-t/τ)) / E.

The expression matches relaxation_time, but the two are not interchangeable: one describes a decay of stress, the other an approach to a bounded strain.

Named chains

MeanFieldHomogenization.Viscoelasticity.zener_maxwell Function

zener_maxwell(E_inf, E_1, tau_1) -> PronyRelaxation

The standard linear solid in its Maxwell (relaxation) representation: an equilibrium spring E_inf in parallel with one Maxwell branch.

zener_kelvin is the same material written the other way round, and maxwell_to_kelvin converts between them exactly — with, for one branch, the classical closed form

MeanFieldHomogenization.Viscoelasticity.zener_kelvin Function

zener_kelvin(E_glassy, E_delayed, tau_1) -> PronyCreep

The standard linear solid in its Kelvin (creep) representation: a spring of modulus E_glassy in series with one Kelvin cell of modulus E_delayed and retardation time tau_1.

See zener_maxwell and kelvin_to_maxwell.

MeanFieldHomogenization.Viscoelasticity.burgers Function

burgers(k_s, eta_s, k_p, eta_p) -> PronyCreep

The Burgers model: a Maxwell unit (k_s, eta_s) in series with a Kelvin cell (k_p, eta_p).

A fluid, so kelvin_to_maxwell returns a two-branch PronyRelaxation with no equilibrium spring — one branch more than the Kelvin form has, which is the degree count of the rational transform. Its relaxation function then matches the classical closed-form cosh/sinh expression to machine precision, which is how the fluid path of the conversion is validated (test/Viscoelasticity/test_prony.jl).

Parameter names follow the ECHOES reference tests/python/creep/fluage_echoes_ijss2013_jsanahuja_relaxBurgers.py.

Fractional elements

MeanFieldHomogenization.Viscoelasticity.ScottBlair Type

ScottBlair(V, α)

The springpot (fractional dashpot, Scott-Blair element), interpolating continuously between a spring (α = 0) and a dashpot (α = 1):

0 < α < 1. The exact pair     is what makes the power-law terms of HuetSayegh and Model2S2P1D analytic in both domains.

The transform has a branch cut along the negative real axis, which the Talbot contours enclose rather than cross — FixedTalbot inverts it at 1e-12, no special handling needed.

MeanFieldHomogenization.Viscoelasticity.FractionalMaxwell Type

FractionalMaxwell(V_a, α, V_b, β)

Two springpots in series:  .

By convention α > β, so the α element is the stiffer one at short times. A fluid.

MeanFieldHomogenization.Viscoelasticity.FractionalKelvin Type

FractionalKelvin(V_a, α, V_b, β)

Two springpots in parallel:   , with α > β.

With β = 0 the slow element is a spring and the model is a solid of equilibrium modulus V_b.

MeanFieldHomogenization.Viscoelasticity.FractionalZener Type

FractionalZener(E_inf, E_0, tau, α)

The fractional standard solid, also known as the Cole-Cole model:

with the one-parameter Mittag-Leffler function. α = 1 recovers zener_maxwell; smaller α broadens the relaxation spectrum, which is what makes the model fit polymers and bitumen where a single exponential cannot.

The time-domain form needs MittagLeffler.jl (a weak dependency); without it, relaxation falls back to numerically inverting the closed-form transform above, which costs a little accuracy and nothing else.

MeanFieldHomogenization.Viscoelasticity.Rabotnov Type

Rabotnov(mu_0, lambda_0, α, β)

The Rabotnov fractional-exponential kernel, in the form ECHOES uses as its analytical benchmark for the ageing pipeline (tests/python/creep/fluage_echoes_maxwell_papier_rabotnov.py):

The Carson transform is elementary

Since    , the Carson transform of that kernel collapses to

with no Mittag-Leffler function anywhere. The model is therefore complete without the MittagLeffler extension: the transform is closed form and the inversion supplies the time domain to about 1e-10. With the extension loaded, relaxation uses the closed form instead — the two agree, which is one of the cross-checks in test/Viscoelasticity/test_rheology.jl.

mu_0 is the glassy modulus (the term vanishes as p → ∞) and   the equilibrium one. A passive material therefore has  : the benchmark of [59] §5 uses μ₀ = 1.7, λ₀ = -0.495, α = -0.46, β = 0.98, and α ∈ (-1, 0) is the usual range — it is α + 1 that must be positive, not α.

Bituminous binders and long-term creep

MeanFieldHomogenization.Viscoelasticity.HuetSayegh Type

HuetSayegh(E00, E0, δ, τ, k, h)

The Huet-Sayegh model: two parabolic (springpot) elements and a spring in series, the whole in parallel with a spring E00.

E00 is the static modulus, E0 the glassy one, and 0 < k < h < 1. A solid: unlike Model2S2P1D there is no series dashpot, so the strain stays bounded.

MeanFieldHomogenization.Viscoelasticity.Model2S2P1D Type

Model2S2P1D(E00, E0, δ, τ_E, k, h, β)

The 2S2P1D model of Di Benedetto and Olard — 2 Springs, 2 Parabolic elements, 1 Dashpot — the reference rheological model for bituminous binders and mixtures:

E00 is the static modulus, E0 the glassy one, 0 < k < h < 1, and the series dashpot β makes the model a fluid.

One model, an exact pair in both domains

Because    , the denominator above is the Laplace-Carson transform of an ordinary function of time,

term by term — see creep_kernel. So 2S2P1D is simultaneously an exact Laplace-Carson model and an exact Volterra model: the same material can be pushed through homogenize_lc and through homogenize_alv with no approximation on the input side, and the two must agree. That is what makes it the headline cross-validation of the two routes.

creep_kernel_law returns packaged as a ViscoLaw for the ageing pipeline.

The field names in older code are swapped

docs/src/applications/bituminous.md used to define its own struct whose fields were named E0 for the static modulus and Einf for the glassy one — the reverse of the convention here and in the literature. The formula and the positional order were already right, so migrating that page is a pure rename with no numerical change; the twelve calibrated parameter sets carry over argument for argument.

MeanFieldHomogenization.Viscoelasticity.creep_kernel Function

creep_kernel(m::Model2S2P1D, t)

The dimensionless creep kernel in the time domain,

the exact term-by-term counterpart of carson_creep_kernel.

MeanFieldHomogenization.Viscoelasticity.carson_creep_kernel Function

carson_creep_kernel(m::Model2S2P1D, p)

The dimensionless denominator of the 2S2P1D transform — the Laplace-Carson transform of creep_kernel.

MeanFieldHomogenization.Viscoelasticity.creep_kernel_law Function

creep_kernel_law(m::Model2S2P1D) -> ViscoLaw

The dimensionless 2S2P1D creep kernel packaged as a scalar :creep ViscoLaw (t, t') ↦ φ(t - t'), ready for the ageing pipeline.

This is the entry point that lets one 2S2P1D object drive both routes. The relaxation modulus in the time domain is

where is the Volterra inverse, not the pointwise one:

m   = Model2S2P1D(1e-7, 1000.0, 2.2, 1.945e-3, 0.22, 0.63, 50.0)
Φ   = trapezoidal_matrix(creep_kernel_law(m), times)     # n × n
R   = m.E00 * I + (m.E0 - m.E00) * volterra_inverse(Φ; block_size = 1)

and that discrete R must agree with relaxation(m, t) obtained by inverting carson_relaxation — the two-route consistency check of test/Viscoelasticity/test_rheology.jl.

Reproduces VM(...) of the ECHOES reference tests/python/creep/modele2S2P1D.py.

MeanFieldHomogenization.Viscoelasticity.LogarithmicCreep Type

LogarithmicCreep(E, C, τ)

The logarithmic creep law used for concrete at long times:

with the exponential integral. A fluid in the sense that the strain grows without bound, though only logarithmically.

The transform is evaluated through SpecialFunctions.expintx, the scaled form , so it stays finite for large pτ where exp(pτ) alone would overflow.

This is the non-ageing skeleton of the ageing law logcompliance in the ECHOES ageing_visco_mat.py; the ageing version, where E, C and τ depend on the loading age, belongs to the time-domain route instead.

Lifting to tensors, and homogenizing

MeanFieldHomogenization.Viscoelasticity.AbstractIsoPairing Type

AbstractIsoPairing

Which pair of scalar quantities an IsoRheology is built from. A trait parameterizing the type rather than a Symbol field, so the choice is resolved at compile time and stays extensible.

Concrete pairings: BulkShear, YoungPoisson.

MeanFieldHomogenization.Viscoelasticity.BulkShear Type

BulkShear

Pairing by bulk and shear: the two channels of the isotropic tensor relax independently, which is the physically natural description and the one every scheme works in internally.

MeanFieldHomogenization.Viscoelasticity.YoungPoisson Type

YoungPoisson

Pairing by Young's modulus and Poisson's ratio, the pair a laboratory reports.

A viscoelastic Poisson ratio is a Laplace-Carson object

The conversion

is a pointwise algebraic identity in the Carson domain only. In the time domain the same relation is a Volterra quotient: ν(t) does not divide, it deconvolves. (That is exactly what the ECHOES reference does with V3k = VE · inv(I - 2Vν) on discretized operators — tests/python/creep/modele2S2P1D.py.)

So a time-dependent Poisson ratio belongs on this side of the fence. A constant nu is unambiguous and is handled as a special case, including for the closed-form time-domain values.

MeanFieldHomogenization.Viscoelasticity.IsoRheology Type

IsoRheology{P, A, B}(a, b)

An isotropic fourth-order viscoelastic model built from two scalar channels a and b, combined according to the pairing P.

Build one with iso_rheology or iso_rheology_E_nu rather than by calling the constructor.

It answers the same generics as a scalar model — carson_relaxation, carson_creep, relaxation, creep, complex_modulus — with TensISO{4,3} values, and converts to a ViscoLaw for the ageing pipeline.

MeanFieldHomogenization.Viscoelasticity.iso_rheology Function

iso_rheology(k_model, mu_model) -> IsoRheology{BulkShear}

Isotropic model whose bulk channel follows k_model and shear channel mu_model, each any AbstractRheology.

A Real in either slot is promoted to a Spring, so the common case of an elastic bulk modulus and a relaxing shear one reads

iso_rheology(2500.0, Model2S2P1D(1e-7, 1000.0, 2.2, 1.945e-3, 0.22, 0.63, 50.0))

Examples

julia> using MeanFieldHomogenization, TensND

julia> m = iso_rheology(zener_maxwell(30.0, 10.0, 1.0), zener_maxwell(10.0, 5.0, 0.5));

julia> C = relaxation(m, 0.0);      # the glassy stiffness

julia> all(isapprox.(get_data(C), (3 * 40.0, 2 * 15.0)))
true

MeanFieldHomogenization.Viscoelasticity.iso_rheology_E_nu Function

iso_rheology_E_nu(E_model, nu) -> IsoRheology{YoungPoisson}

Isotropic model given by a Young's-modulus channel and a Poisson ratio, where nu is either a constant Real or itself an AbstractRheology.

See YoungPoisson for why a time-dependent nu only makes unambiguous sense in the Laplace-Carson domain.

MeanFieldHomogenization.Viscoelasticity.homogenize_lc Function

homogenize_lc(build_cell, scheme, property = :C; p)
homogenize_lc(build_cell, scheme, property = :C; times, method = DEFAULT_INVERSION, kw...)

Homogenize a non-ageing viscoelastic composite through the correspondence principle, and — in the second form — bring the answer back to the time domain.

build_cell(p) is a closure returning the homogenization cell with every phase property evaluated at the Carson variable p. Any cell works: an RVE, a LayeredSphere, a Laminate, a nested chain of them. Extra keyword arguments are forwarded to homogenize.

Example

matrix = iso_rheology(Spring(2500.0), Model2S2P1D(1e-7, 1000.0, 2.2, 1.945e-3, 0.22, 0.63, 50.0))
aggreg = iso_rheology(Spring(30000.0), Spring(22000.0))

function cell(p)
    rve = RVE()
    add_phase!(rve, :MASTIC, Ellipsoid(1.0), Dict(:C => carson_relaxation(matrix, p)); fraction = :rest)
    add_phase!(rve, :AGG, Ellipsoid(1.0), Dict(:C => carson_relaxation(aggreg, p));
               fraction = 0.35)
    return rve
end

E_star = homogenize_lc(cell, MoriTanaka(), :C; p = im * 2π * 10)     # 10 Hz
R_of_t = homogenize_lc(cell, MoriTanaka(), :C; times = exp10.(-4:0.1:2))

Choosing the inversion method

Each time costs one inversion, and each inversion costs N full homogenizations, so the choice is a real one:

• FixedTalbot (the default) — most accurate, 24 complex homogenizations per time;

• GaverStehfest — 16 real homogenizations per time. Real arithmetic is roughly twice as cheap per evaluation, and it is the only way to use SelfConsistent(algorithm = NewtonDefault()) here, whose ForwardDiff Jacobian cannot carry a Dual over a complex scalar. The price is about five significant digits instead of twelve;

• DeHoog — one node set serves a whole block of times, so a multi-decade grid costs far fewer homogenizations than any per-point method. This is usually the right answer for a master curve.

Relation to the time-domain route

For a non-ageing material this and homogenize_alv compute the same thing by disjoint means — no shared code — so agreement between them is a real check on both. ViscoLaw(model) turns the same model objects into the kernels homogenize_alv needs.

Optional Mittag-Leffler support

MeanFieldHomogenization.Viscoelasticity._mittag_leffler Function

_mittag_leffler(a, b, z) -> value or `nothing`

The two-parameter Mittag-Leffler function .

Returns nothing when MittagLeffler.jl is not loaded, which is the signal callers use to fall back on inverting the closed-form Carson transform:

ml = _mittag_leffler(α, one(α), -(t / τ)^α)
ml === nothing && return inverse_carson(p -> carson_relaxation(m, p), t, method)

Load the package to switch the closed forms on:

using MittagLeffler        # activates MeanFieldHomogenizationMittagLefflerExt

The extension declines more than it accepts

MittagLeffler.jl v1.0.0 silently returns 1.0 for 1 < a < 2 and small |z| — mittleff(1.3, 1.0, -0.5) gives 1.0 where the defining series gives 0.6330079. The extension therefore answers only for 0 < a ≤ 1 and b == 1, the domain checked against that series.

FractionalZener (0 < α ≤ 1) and Rabotnov with the physical α ∈ (-1, 0) — whose Mittag-Leffler order α + 1 then lies in (0, 1) — are both inside that domain and get their closed forms. A Rabotnov kernel with α > 0 is not, and falls back on the inversion, which agrees with the closed form to about 1e-10 wherever both can be evaluated.

This replaces a PyCall detour

Before this extension, the Rabotnov benchmark reached E_{a,b} through PyCall and an external Python module living outside the repository (scripts/52_rabotnov_mittag_leffler.jl, and §11 of the viscoelasticity manual). That is no longer needed: the transform is closed-form and the time-domain reference is a registered Julia package away.