Fluidity-Saramito EVP Model

Quick Reference

  • Use when: Elastoviscoplastic materials with thixotropy

  • Parameters: 10 (minimal) or 12 (full coupling): \(G\), \(\eta_s\), \(\tau_{y0}\), \(K_{HB}\), \(n_{HB}\), \(f_{\text{age}}\), \(f_{\text{flow}}\), \(t_a\), \(b\), \(n_{\text{rej}}\); full adds \(\tau_{y,\text{coupling}}\), \(m_{\text{yield}}\)

  • Key equation: \(\boldsymbol{\tau} + \lambda \stackrel{\nabla}{\boldsymbol{\tau}} = 2\eta_p \mathbf{D}\)

  • Test modes: flow_curve, startup, creep, oscillation, LAOS

  • Material examples: Carbopol, hair gel, mayonnaise

Model Summary

Model Class

FluiditySaramitoLocal, FluiditySaramitoNonlocal

Physics

Elastoviscoplastic with thixotropic fluidity

Coupling Modes

"minimal", "full"

Protocols (Local)

FLOW_CURVE, CREEP, RELAXATION, STARTUP, OSCILLATION, LAOS

Protocols (Nonlocal)

FLOW_CURVE, CREEP, STARTUP

Key Features

Tensorial stress, Von Mises yield, normal stresses, shear banding

Import:

from rheojax.models.fluidity import FluiditySaramitoLocal, FluiditySaramitoNonlocal

Basic Usage:

# Minimal coupling (simplest, most identifiable)
model = FluiditySaramitoLocal(coupling="minimal")
model.fit(gamma_dot, sigma, test_mode='flow_curve')

# Full coupling (aging yield stress)
model = FluiditySaramitoLocal(coupling="full")
model.fit(gamma_dot, sigma, test_mode='flow_curve')

Notation Guide

Symbol

Description

Units

\(\boldsymbol{\tau}\)

Deviatoric stress tensor

Pa

\(|\boldsymbol{\tau}|\)

Von Mises equivalent stress

Pa

\(\tau_y\)

Yield stress

Pa

\(f\)

Fluidity

1/(Pa·s)

\(\lambda\)

Relaxation time (= 1/f)

s

\(G\)

Elastic modulus

Pa

\(\dot{\gamma}\)

Shear rate

1/s

\(\alpha\)

Plasticity parameter

dimensionless

\(\xi\)

Cooperativity length

m

Overview

The Fluidity-Saramito Elastoviscoplastic (EVP) model combines three key physical mechanisms:

  1. Viscoelasticity: Upper-convected Maxwell framework with elastic recoil, storage modulus \(G'\), and first normal stress difference \(N_1\).

  2. Viscoplasticity: True Von Mises yield surface with Herschel-Bulkley plastic flow above yield.

  3. Thixotropy: Time-dependent aging (structural build-up at rest) and shear rejuvenation (flow-induced breakdown) via fluidity evolution.

The model captures complex behaviors including:

  • Stress overshoot in startup that increases with waiting time

  • Creep bifurcation at the yield stress (bounded vs unbounded flow)

  • Non-exponential stress relaxation

  • Shear banding in spatially-resolved (nonlocal) variant

Physical Foundations

Upper-Convected Maxwell Framework

The stress evolution follows the upper-convected Maxwell model with plasticity:

\[\lambda \overset{\nabla}{\boldsymbol{\tau}} + \alpha(\boldsymbol{\tau})\boldsymbol{\tau} = 2\eta_p \mathbf{D}\]

where:

  • \(\lambda = 1/f\) is the fluidity-dependent relaxation time

  • \(\overset{\nabla}{\boldsymbol{\tau}}\) is the upper-convected derivative

  • \(\alpha = \max(0, 1 - \tau_y/|\boldsymbol{\tau}|)\) is the Von Mises plasticity

  • \(\eta_p = G/f\) is the polymeric viscosity

  • \(\mathbf{D}\) is the rate of deformation tensor

Von Mises Yield Criterion

The plasticity parameter \(\alpha\) activates plastic flow only when the Von Mises equivalent stress exceeds the yield stress:

\[\alpha = \max\left(0, 1 - \frac{\tau_y}{|\boldsymbol{\tau}|}\right)\]

where the Von Mises stress is:

\[|\boldsymbol{\tau}| = \sqrt{\frac{1}{2}\boldsymbol{\tau}:\boldsymbol{\tau}}\]

Fluidity Evolution

The fluidity evolves via competing aging and rejuvenation:

\[\frac{df}{dt} = \frac{f_\text{age} - f}{t_a} + b|\dot{\gamma}|^{n_\text{rej}}(f_\text{flow} - f)\]

where:

  • \(f_\text{age}\): Equilibrium fluidity at rest (aged state)

  • \(f_\text{flow}\): High-shear fluidity limit (rejuvenated state)

  • \(t_a\): Aging timescale

  • \(b\): Rejuvenation amplitude

  • \(n_\text{rej}\): Rejuvenation rate exponent

Coupling Modes

Minimal Coupling (coupling="minimal"):

  • Relaxation time: \(\lambda = 1/f\)

  • Yield stress: \(\tau_y = \tau_{y0}\) (constant)

  • Fewer parameters, easier to identify

Full Coupling (coupling="full"):

  • Relaxation time: \(\lambda = 1/f\)

  • Yield stress: \(\tau_y(f) = \tau_{y0} + a_y/f^m\)

  • Captures aging yield stress (stronger when aged)

Coupling Architecture Design

The Fluidity-Saramito model offers systematic control over how fluidity couples to mechanical properties. Understanding these coupling choices is essential for capturing specific physical behaviors.

Three Coupling Architectures:

  1. Minimal Coupling (coupling="minimal"):

    • Only relaxation time depends on fluidity: \(\lambda = 1/f\)

    • Yield stress constant: \(\tau_y = \tau_{y0}\)

    • Fewest parameters, most identifiable from data

    • Use when: Standard thixotropic EVP, no aging-dependent yield

  2. Aging Yield Coupling (coupling="full"):

    • Relaxation time: \(\lambda = 1/f\)

    • Yield stress increases with aging (lower fluidity):

      \[\tau_y(f) = \tau_{y,\min} + \Delta\tau_y \left(\frac{f_*}{f + f_*}\right)^m\]

    • Captures wait-time dependent yield stress

    • Use when: Materials that strengthen significantly at rest (waxy crude oils, cement, greases)

  3. Dissipation-Consistent Driving (advanced, not yet implemented):

    • Fluidity evolution driven by plastic dissipation \(|\boldsymbol{\tau}:\mathbf{D}|\) instead of kinematic shear rate \(|\dot{\gamma}|\):

      \[\frac{\partial f}{\partial t} = \frac{f_{\rm eq} - f}{\tau_{\rm age}} + b \left(\frac{|\boldsymbol{\tau}:\mathbf{D}|}{\tau_*}\right)^n (f_{\rm flow} - f)\]

    • Thermodynamically consistent formulation

    • Naturally bounds rejuvenation rate during elastic loading

    • Use when: Strong elastic effects, modeling viscoelastic recoil

Coupling Design Checklist:

When selecting coupling options, consider these “knobs” in order of priority:

Knob

Parameter

Effect

Recommendation

1. \(\lambda(f)\)

\(\lambda = 1/f\)

Thixotropic viscosity, stress relaxation rate

Always recommended — core physics

2. \(\tau_y(f)\)

Aging yield coupling

Wait-time dependent yield stress

Use if yield stress varies with rest time

3. G(f)

Elastic modulus coupling

Structural-dependent stiffness

Use sparingly — hard to identify, often negligible

Driving Term Options:

Driving Term

Equation

Physical Interpretation

Kinematic (default)

\(a|\dot{\gamma}|^n(f_{\rm flow} - f)\)

Rejuvenation from deformation rate

Energetic/Dissipative

\(b(|\boldsymbol{\tau}:\mathbf{D}|/\tau_*)^n(f_{\rm flow} - f)\)

Rejuvenation from plastic work

Plastic-only

\(b\alpha|\boldsymbol{\tau}:\mathbf{D}|^n(f_{\rm flow} - f)\)

Only above yield (\(\alpha > 0\))

What the Coupled Model Uniquely Predicts:

Compared to plain Saramito (no fluidity):

  • History-dependent flow curves: Same \(\dot{\gamma}\) gives different \(\sigma\) depending on shear history

  • Waiting-time dependence: Stress overshoot and yield stress increase with rest time

  • Non-monotonic startup: Peak stress before steady state

  • Creep bifurcation: Bounded vs unbounded flow at the yield stress

Compared to plain fluidity (scalar stress):

  • Tensorial stress state: Full \(\boldsymbol{\tau}_{ij}\) evolution

  • Normal stress difference \(N_1\): Rod climbing, die swell

  • Realistic elastic loading: Proper rate-dependent elastic response

  • Von Mises yield criterion: True 3D yielding behavior

Stress-Driven Fluidity Evolution

An alternative to kinematic (shear-rate) driving is stress-driven rejuvenation, where fluidity evolution depends on the deviatoric stress magnitude rather than the deformation rate.

Stress-Driven Evolution Equation:

\[\frac{df}{dt} = \frac{f_{\rm age} - f}{t_a} + a \left(\frac{|\boldsymbol{\tau}_d|}{\tau_y}\right)^{n_{\rm rej}} (f_{\rm flow} - f)\]

where \(|\boldsymbol{\tau}_d| = \sqrt{\frac{1}{2}\boldsymbol{\tau}_d:\boldsymbol{\tau}_d}\) is the deviatoric stress magnitude.

Feedback Loop Dynamics:

Stress-driven rejuvenation creates a self-regulating feedback loop:

  1. Elastic loading phase: Stress builds up, \(|\boldsymbol{\tau}_d|\) increases, rejuvenation term activates

  2. Fluidity increase: Higher \(f\) means faster stress relaxation (\(\alpha f \boldsymbol{\tau}\) term)

  3. Stress saturation: Relaxation limits further stress growth, which in turn limits rejuvenation

  4. Self-limiting equilibrium: System naturally finds steady state

This contrasts with kinematic driving where \(|\dot{\gamma}|\) is externally imposed and can lead to unbounded rejuvenation during rapid startup.

When to Use Stress-Driven vs Kinematic:

Driving Type

Use When

Caution

Kinematic \(|\dot{\gamma}|\)

Standard rheometry, rate-controlled tests

May over-predict rejuvenation during elastic regime

Stress \(|\boldsymbol{\tau}_d|/\tau_y\)

Strong elastic effects, stress-controlled tests

Requires stress evolution equations

Plastic work \(|\boldsymbol{\tau}:\mathbf{D}|\)

Thermodynamic consistency required

Most complex implementation

Physical Interpretation:

Stress-driven rejuvenation captures the physics that microstructural breakdown occurs when the material is under stress, not merely when it deforms. This is particularly relevant for:

  • Materials with significant elastic strain before yielding

  • Stress-controlled creep tests

  • Understanding the viscosity bifurcation near yield stress

  • Materials where elastic recoil is important

Temporal Evolution Regimes

Understanding the transient response of the Fluidity-Saramito model requires tracking the coupled evolution of stress tensor components and fluidity through distinct temporal regimes.

Startup Flow Timeline:

During startup shear from rest (at constant \(\dot{\gamma}\)), the material passes through characteristic stages:

  1. Initial state (\(t = 0\)):

    • Fluidity at aged value: \(f = f_{\rm aged}\) (low, from prior rest)

    • Relaxation time large: \(\lambda = 1/f \gg 1\)

    • Stress tensor: \(\boldsymbol{\tau} = \mathbf{0}\)

  2. Elastic loading (\(t < t_{\rm yield}\), typically \(\gamma < \gamma_y \sim 0.01-0.1\)):

    • Stress builds approximately as: \(\tau_{xy} \approx G \cdot \gamma = G \dot{\gamma} t\)

    • Von Mises stress: \(|\boldsymbol{\tau}| \approx \tau_{xy}\) (shear-dominated)

    • Plasticity inactive: \(\alpha = 0\) (below yield)

    • Fluidity unchanged (no flow → no rejuvenation): \(df/dt \approx (f_{\rm age} - f)/t_a \approx 0\)

  3. Yield onset (\(|\boldsymbol{\tau}| \to \tau_y\)):

    • Plasticity activates: \(\alpha = \max(0, 1 - \tau_y/|\boldsymbol{\tau}|) > 0\)

    • Stress growth rate slows (relaxation now competes with loading)

    • For full coupling: \(\tau_y(f_{\rm aged})\) may exceed \(\tau_{y0}\) significantly

  4. Peak overshoot (\(t = t_{\rm peak}\), typically \(\gamma \sim 0.1-1\)):

    • Maximum stress: \(\sigma_{\rm max} = \tau_y(f_{\rm aged}) + G \gamma_y + \text{viscous contribution}\)

    • Balance point: elastic loading rate \(\approx\) plastic dissipation rate

    • Fluidity beginning to increase (rejuvenation activating)

  5. Rejuvenation-dominated decay (\(t > t_{\rm peak}\)):

    • Fluidity increases: \(df/dt = b|\dot{\gamma}|^{n_{\rm rej}}(f_{\rm flow} - f) > 0\)

    • Relaxation time decreases: \(\lambda = 1/f \downarrow\)

    • Stress decays toward steady state

    • For full coupling: \(\tau_y(f)\) also decreases

  6. Steady state (\(t \gg t_a, t_{\rm peak}\)):

    • Aging = rejuvenation: \((f_{\rm age} - f)/t_a = b|\dot{\gamma}|^{n_{\rm rej}}(f - f_{\rm flow})\)

    • Stress: \(\sigma_{ss} = \tau_y(f_{ss}) + K_{\rm HB}\dot{\gamma}^{n_{\rm HB}}\)

Creep Bifurcation Dynamics:

Under constant applied stress \(\sigma_{\rm applied}\):

  1. Sub-yield (\(\sigma_{\rm applied} < \tau_y\)):

    • Initial elastic strain: \(\gamma_e = \sigma_{\rm applied}/G\)

    • Plasticity inactive: \(\alpha = 0\)

    • No rejuvenation (no flow): fluidity decreases (aging)

    • \(f \to f_{\rm age}\), \(\tau_y(f) \to \tau_y(f_{\rm age}) > \tau_y(f_{\rm flow})\)

    • Material stiffens: bounded strain, solid-like arrest

  2. Above yield (\(\sigma_{\rm applied} > \tau_y\)):

    • Initial elastic strain plus viscoplastic flow

    • Plasticity active: \(\alpha > 0 \Rightarrow\) positive \(\dot{\gamma}\)

    • Positive feedback loop:

      • Flow → rejuvenation → \(f \uparrow\)

      • \(f \uparrow \Rightarrow \lambda \downarrow \Rightarrow\) faster flow

      • For full coupling: \(f \uparrow \Rightarrow \tau_y \downarrow \Rightarrow\) easier yielding

    • Result: avalanche-like transition to steady flow

  3. Near yield (\(\sigma_{\rm applied} \approx \tau_y\)):

    • Delayed yielding: Long induction time before flow accelerates

    • Metastable: small perturbations determine bounded vs unbounded outcome

    • Duration of plateau: \(t_{\rm delay} \sim t_a \cdot \ln(\Delta f / \epsilon)\) where \(\epsilon\) is initial departure from equilibrium

Flow Cessation (Relaxation) Dynamics:

When shear stops (\(\dot{\gamma} \to 0\)) from steady flow:

  1. Immediate response (\(t = 0^+\)):

    • Rejuvenation term vanishes: \(b|\dot{\gamma}|^n(f_{\rm flow} - f) \to 0\)

    • Stress locked in current configuration

    • Elastic recoil (partial strain recovery) on timescale \(\sim \lambda\)

  2. Aging phase (\(t > 0\)):

    • Pure aging: \(df/dt = (f_{\rm age} - f)/t_a < 0\)

    • Fluidity decays: \(f(t) = f_{\rm age} + (f_{ss} - f_{\rm age})e^{-t/t_a}\)

    • Relaxation time grows: \(\lambda(t) = 1/f(t) \uparrow\)

  3. Stress relaxation:

    • Below yield (\(|\boldsymbol{\tau}| < \tau_y\)): Stress frozen (elastic solid)

    • At/above yield: Slow viscoplastic relaxation on growing timescale

    • Non-exponential: \(\tau(t) \sim |\boldsymbol{\tau}(0)| \cdot \exp\left(-\int_0^t \alpha(s) f(s) ds\right)\)

  4. Long-time behavior:

    • \(f \to f_{\rm age}\) exponentially with time constant \(t_a\)

    • For full coupling: \(\tau_y(f) \to \tau_y(f_{\rm age})\) (strengthening)

    • Residual stress relaxes slowly if \(|\boldsymbol{\tau}| > \tau_y\)

Regime Summary Diagram:

Startup:
───────────────────────────────────────────────────────────────►
│     Elastic    │   Yield   │    Peak     │  Rejuvenation  │ SS
│     loading    │   onset   │   stress    │     decay      │
f: ▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔───────────────────────────── → f_ss
σ:              ╱╲
              ╱    ╲─────────────────────────────────────────── → σ_ss
            ╱
         ╱
▁▁▁▁▁▁╱
0    γ_y                                                    t →

Mathematical Formulation

Governing Equations

Component Equations (Simple Shear)

For simple shear with velocity gradient \(\mathbf{L} = \dot{\gamma}\mathbf{e}_x\mathbf{e}_y\):

\[\begin{split}\frac{d\tau_{xx}}{dt} &= 2\dot{\gamma}\tau_{xy} - \alpha f \tau_{xx} \\ \frac{d\tau_{yy}}{dt} &= -\alpha f \tau_{yy} \\ \frac{d\tau_{xy}}{dt} &= \dot{\gamma}\tau_{yy} + G\dot{\gamma} - \alpha f \tau_{xy}\end{split}\]

The first normal stress difference is:

\[N_1 = \tau_{xx} - \tau_{yy}\]

At steady state in simple shear, this scales as \(N_1 \sim \lambda \dot{\gamma} \tau_{xy}\).

Steady-State Flow Curve

At steady state, the model reduces to Herschel-Bulkley form:

\[\sigma = \tau_y + K_\text{HB}\dot{\gamma}^{n_\text{HB}}\]

with fluidity-dependent parameters when using full coupling.

Validity and Assumptions

Model Assumptions

  1. Tensorial stress description: Uses full stress tensor \(\boldsymbol{\tau}\), suitable for 3D flows

  2. Von Mises yield criterion: Isotropic yielding based on stress magnitude \(|\boldsymbol{\tau}|\)

  3. Upper-convected Maxwell framework: Appropriate for entangled polymers and structured fluids

  4. Thixotropic fluidity evolution: Single scalar internal variable \(f\) tracks structure

  5. Isothermal: Temperature effects not explicitly modeled

Applicability

The model is most appropriate for:

  • Elastoviscoplastic fluids: Materials exhibiting elastic recoil, yield stress, and time-dependent viscosity

  • Concentrated emulsions: Mayonnaise, creams, cosmetics

  • Polymer gels: Carbopol, hydrogels with yield stress

  • Pastes and slurries: Cement, drilling muds, waxy crude oils

  • Complex loading: Protocols requiring tensorial stress (extensional flow, normal stresses)

Limitations

No shear banding (local model):

The homogeneous (local) model cannot capture spatial heterogeneity. For shear-banded flows, use the FluiditySaramitoNonlocal variant.

Single structural variable:

Real materials may have multiple structural timescales or multi-mode thixotropy. Consider multi-lambda extensions for complex aging.

Isotropic yielding:

The Von Mises criterion assumes isotropic yield surface. Anisotropic materials may require tensorial yield criteria.

Small strain assumption:

The UCM framework assumes affine deformation. Not suitable for materials with wall slip or non-affine microstructure.

What You Can Learn

From fitting Fluidity-Saramito EVP to experimental data, you can extract insights about elastoviscoplasticity, thixotropic coupling, and normal stress generation in yield-stress materials.

Parameter Interpretation

f (Fluidity):

Time-dependent inverse relaxation time controlling both viscosity \(\eta_p = G/f\) and relaxation time \(\lambda = 1/f\) in tensorial UCM framework. For graduate students: Unlike scalar fluidity models, here f couples to tensorial stress evolution: \(\lambda \cdot \overset{\nabla}{\boldsymbol{\tau}} + \alpha(\boldsymbol{\tau}) \boldsymbol{\tau} = 2\eta_p \mathbf{D}\) where \(\lambda = 1/f\). Evolution: \(df/dt = (f_{\text{age}} - f)/t_a + b|\dot{\gamma}|^{n_{\text{rej}}}(f_{\text{flow}} - f)\). Plasticity \(\alpha = \max(0, 1 - \tau_y/|\boldsymbol{\tau}|)\) activates only when Von Mises stress \(|\boldsymbol{\tau}| = \sqrt{\frac{1}{2}\boldsymbol{\tau}:\boldsymbol{\tau}} > \tau_y\). Normal stresses \(N_1 \sim \lambda \dot{\gamma} \tau_{xy} \sim \dot{\gamma}/f^2\) (aged materials with low f exhibit larger \(N_1\)). For practitioners: \(f_{\text{age}} \approx 10^{-6}\) to \(10^{-3}\) s-1 (solid-like), \(f_{\text{flow}} \approx 10^{-2}\) to 1 s-1 (liquid-like). Measure via startup overshoot magnitude (larger overshoot = lower initial \(f\) after aging).

G (Elastic Modulus):

Elastic stiffness in UCM backbone, controlling both stress buildup and normal stress generation. For graduate students: Sets stress scale and Weissenberg number \(\text{Wi} = \lambda \dot{\gamma} = \dot{\gamma}/f\). Normal stress scaling: \(N_1/\tau_{xy} \sim \text{Wi} \sim \dot{\gamma}/f\). For UCM, \(N_2 = 0\). Startup overshoot: \(\sigma_{\text{peak}} \sim \tau_y + G\gamma_y\) where \(\gamma_y\) is yield strain. For practitioners: Extract from initial slope in startup or from \(G'\) plateau in SAOS. Typical: \(G = 10^2\text{--}10^4\) Pa (soft colloids), \(10^4\text{--}10^6\) Pa (polymer gels).

\(\tau_y\) (Yield Stress):

Von Mises yield criterion threshold. In minimal coupling, constant. In full coupling: \(\tau_y(f) = \tau_{y0} + a_y/f^m\) (aged materials are stronger). For graduate students: Plasticity parameter \(\alpha = \max(0, 1 - \tau_y/|\tau|)\) controls plastic dissipation term. \(\alpha = 0\) below yield (elastic), \(\alpha > 0\) above yield (viscoplastic flow). Full coupling exponent \(m\) typically 0.3-1.0 captures aging-induced hardening. For practitioners: Measure from flow curve low-shear plateau or creep bifurcation stress. Full coupling appropriate if yield stress increases significantly after aging (wait-time dependent startup tests).

t_a, b, n_rej (Fluidity Evolution Parameters):

Control aging (\(t_a\)) and rejuvenation (\(b\), \(n_{\text{rej}}\)) kinetics analogous to local fluidity model. For graduate students: \(df/dt = (f_{\text{age}} - f)/t_a + b|\dot{\gamma}|^{n_{\text{rej}}}(f_{\text{flow}} - f)\). Characteristic shear rate: \(\dot{\gamma}_c \sim (1/(bt_a))^{1/n_{\text{rej}}}\). Stress overshoot magnitude and position depend on \(t_a\) (waiting time scaling), \(b\) and \(n_{\text{rej}}\) (breakdown kinetics). For practitioners: Extract from startup test family at different wait times. Typical: \(t_a = 10\text{--}1000\) s, \(b = 0.1\text{--}10\), \(n_{\text{rej}} = 0.5\text{--}1.5\). Longer \(t_a\) means more pronounced thixotropy.

Material Classification

Material Classification from Fluidity-Saramito Parameters

Parameter Range

Material Behavior

Typical Materials

Processing Implications

\(N_1/\tau_{xy} < 0.1\), minimal coupling

Weakly elastic EVP

Carbopol gels, pastes

Dominant viscoplasticity, minimal normal stresses

\(N_1/\tau_{xy} = 0.1\text{--}1\), minimal coupling

Moderate elasticity

Emulsions, soft colloids, cosmetics

Significant rod climbing, moderate die swell

\(N_1/\tau_{xy} > 1\), minimal coupling

Strongly elastic EVP

Polymer gels, fiber suspensions

Strong Weissenberg effect, edge fracture risk

Full coupling: m > 0.5

Strong aging-yield coupling

Waxy crude oils, cement

Wait-time dependent yield, restart challenges

\(t_a > 100\) s, full coupling

Extreme thixotropy with elasticity

Waxy crude oils, drilling muds

Long-term aging, complex restart dynamics

Parameters

Core Parameters

Parameter

Description

Units

Default

Bounds

G

Elastic modulus

Pa

1e4

[1e1, 1e8]

tau_y0

Base yield stress

Pa

100

[0.1, 1e5]

K_HB

HB consistency index

Pa·s^n

50

[1e-2, 1e5]

n_HB

HB flow exponent

—

0.5

[0.1, 1.5]

Fluidity Parameters

Parameter

Description

Units

Default

Bounds

f_age

Aging fluidity limit

1/(Pa·s)

1e-6

[1e-12, 1e-2]

f_flow

Flow fluidity limit

1/(Pa·s)

1e-2

[1e-6, 1.0]

t_a

Aging timescale

s

10

[0.01, 1e5]

b

Rejuvenation amplitude

—

1.0

[0, 1e3]

n_rej

Rejuvenation exponent

—

1.0

[0.1, 3.0]

Full Coupling Parameters

Only active when coupling="full":

Parameter

Description

Units

Default

Bounds

tau_y_coupling

Yield stress coupling

Pa·(Pa·s)^m

1.0

[0, 1e4]

m_yield

Yield stress exponent

—

0.5

[0.1, 2.0]

Nonlocal Parameters

Only for FluiditySaramitoNonlocal:

Parameter

Description

Units

Default

Bounds

xi

Cooperativity length

m

1e-5

[1e-7, 1e-2]

Parameter Interpretation by Material

Concentrated Emulsions (mayonnaise, cosmetics):

  • \(\tau_y \sim 10-100\) Pa

  • \(n_\text{HB} \sim 0.3-0.5\)

  • \(t_a \sim 10-100\) s

Polymer Gels (carbopol, hydrogels):

  • \(\tau_y \sim 1-50\) Pa

  • \(n_\text{HB} \sim 0.4-0.6\)

  • \(t_a \sim 1-1000\) s (depends on concentration)

Cement/Concrete:

  • \(\tau_y \sim 100-1000\) Pa

  • \(n_\text{HB} \sim 0.2-0.4\)

  • \(t_a \sim 100-10000\) s (hydration-dependent)

Drilling Muds:

  • \(\tau_y \sim 5-50\) Pa

  • \(n_\text{HB} \sim 0.3-0.7\)

  • \(t_a \sim 10-1000\) s

Usage

(Usage examples already present in the file)

See Also

Fitting Guidance

Step-by-Step Example

from rheojax.models.fluidity import FluiditySaramitoLocal
import numpy as np

# 1. Flow curve fitting
model = FluiditySaramitoLocal(coupling="minimal")
model.fit(gamma_dot, sigma, test_mode='flow_curve')

# 2. Startup fitting (refines G and fluidity parameters)
model.fit(t, sigma_startup, test_mode='startup', gamma_dot=1.0)

# 3. Bayesian inference for uncertainty
result = model.fit_bayesian(
    gamma_dot, sigma,
    test_mode='flow_curve',
    num_warmup=1000,
    num_samples=2000,
    num_chains=4,
)

# 4. Get credible intervals
intervals = model.get_credible_intervals(result.posterior_samples)
for param, stats in intervals.items():
    print(f"{param}: {stats['mean']:.2f} [{stats['lower']:.2f}, {stats['upper']:.2f}]")

Troubleshooting

Issue

Solution

No stress overshoot

Increase b (rejuvenation) or decrease t_a (aging time)

Overshoot too large

Decrease b or increase f_age

Flow curve too flat

Decrease n_HB (more shear-thinning)

Poor creep fit

Check tau_y0 against bifurcation point in data

Bayesian divergences

Use NLSQ warm-start, increase num_warmup

Parameter Estimation from Multiple Protocols

For reliable parameter estimation in the Fluidity-Saramito model, a sequential multi-protocol approach is essential. Each protocol isolates different parameters, enabling progressive refinement.

Sequential Fitting Workflow:

Step

Protocol

Parameters Constrained

Notes

1

SAOS (small amplitude oscillation)

\(G\) from \(G'\) plateau; initial \(f\) from \(\tan\delta\)

Linear regime only. \(G = G'(\omega \to \infty)\) or plateau

2

Flow curve (steady shear)

\(\tau_{y0}\), \(K_{\rm HB}\), \(n_{\rm HB}\)

Fit Herschel-Bulkley. Fixes steady-state rheology

3

Startup tests (multiple \(\dot{\gamma}\))

\(t_a\) (from peak timing), \(b\), \(n_{\rm rej}\) (from peak magnitude)

Use waiting time series. Overshoot position → kinetics

4

Creep tests (\(\sigma\) near \(\tau_y\))

Validate \(\tau_{y0}\) from bifurcation point

\(\sigma < \tau_y\): bounded; \(\sigma > \tau_y\): unbounded. Confirms yield

Step 1: SAOS Analysis

Extract elastic modulus \(G\) and estimate initial fluidity:

# From SAOS data (G', G'' vs omega)
G_prime_plateau = np.max(G_prime)  # High-frequency plateau
G = G_prime_plateau

# Estimate fluidity from loss tangent at intermediate frequency
# tan(δ) = G''/G' ≈ 1/(ωλ) = f·ω for Maxwell-like response
omega_mid = omega[len(omega)//2]
tan_delta_mid = G_double_prime[len(omega)//2] / G_prime[len(omega)//2]
f_initial_estimate = tan_delta_mid * omega_mid

print(f"Estimated G: {G:.1f} Pa")
print(f"Estimated initial f: {f_initial_estimate:.2e} 1/(Pa·s)")

Step 2: Flow Curve Fitting

Fix HB steady-state parameters:

from rheojax.models.fluidity import FluiditySaramitoLocal

model = FluiditySaramitoLocal(coupling="minimal")

# Fix G from SAOS, fit HB parameters only
model.parameters.set_value("G", G, fixed=True)
model.fit(gamma_dot, sigma, test_mode='flow_curve')

tau_y0 = model.parameters.get_value("tau_y0")
K_HB = model.parameters.get_value("K_HB")
n_HB = model.parameters.get_value("n_HB")

Step 3: Startup Test Analysis

Extract kinetic parameters from stress overshoot:

# Startup at multiple waiting times
wait_times = [1, 10, 100, 1000]  # seconds
gamma_dot_startup = 1.0

overshoot_data = []
for t_wait in wait_times:
    _, stress, _ = model.simulate_startup(t, gamma_dot_startup, t_wait=t_wait)
    sigma_max = np.max(stress)
    t_peak = t[np.argmax(stress)]
    sigma_ss = stress[-1]
    overshoot_data.append({
        't_wait': t_wait,
        'sigma_max': sigma_max,
        't_peak': t_peak,
        'overshoot_ratio': sigma_max / sigma_ss
    })

# t_a controls: log(sigma_max) vs log(t_wait) slope
# b controls: overall overshoot magnitude
# n_rej controls: dependence on shear rate

Parameter Sensitivity Guide:

Parameter

Most Sensitive Observables

Identifiability Notes

\(G\)

Initial slope in startup; \(G'\) plateau in SAOS

Well-identified from SAOS or early startup

\(\tau_{y0}\)

Flow curve low-shear limit; creep bifurcation

Requires low \(\dot{\gamma}\) data (< 0.01 s-1)

\(K_{\rm HB}\), \(n_{\rm HB}\)

Flow curve shape

Often correlated; fix one if uncertain

\(t_a\)

Overshoot position vs wait time; relaxation timescale

Best from wait-time series

\(b\)

Overshoot magnitude; rejuvenation rate

Correlated with \(n_{\rm rej}\). Fix one first

\(n_{\rm rej}\)

Shear-rate dependence of rejuvenation

Often fixed at 1.0 initially

\(f_{\rm age}\), \(f_{\rm flow}\)

Aged vs flowing viscosity limits

Use extreme conditions (long rest, high shear)

Spatially-Resolved Data (Nonlocal Variant):

For the nonlocal model, velocity profiles provide direct access to \(\xi\):

  1. Velocity profile measurement: Ultrasound velocimetry (USV), MRI, or PIV

  2. Shear rate extraction: \(\dot{\gamma}(y) = dv/dy\)

  3. Fluidity profile: \(f(y) = \dot{\gamma}(y)/\sigma\) (uniform stress)

  4. \(\xi\) extraction: Fit boundary layer decay of \(f(y)\) from walls

# From velocity profile data: v(y) and known stress Sigma
gamma_dot_profile = np.gradient(v_profile, y_grid)
f_profile = gamma_dot_profile / Sigma

# Fit exponential boundary layer: f(y) = f_bulk + (f_wall - f_bulk)*exp(-y/xi)
from scipy.optimize import curve_fit

def boundary_layer(y, f_bulk, f_wall, xi):
    return f_bulk + (f_wall - f_bulk) * np.exp(-y / xi)

# Fit near wall (y < H/4)
mask = y_grid < H / 4
popt, _ = curve_fit(boundary_layer, y_grid[mask], f_profile[mask],
                    p0=[f_profile[-1], f_profile[0], 50e-6])
xi_extracted = popt[2]

Usage Examples

Basic Flow Curve Fitting

from rheojax.models.fluidity import FluiditySaramitoLocal
import numpy as np

# Generate synthetic data
gamma_dot = np.logspace(-2, 2, 30)
sigma_data = 100 + 50 * gamma_dot**0.5  # HB-like

# Fit model
model = FluiditySaramitoLocal(coupling="minimal")
model.fit(gamma_dot, sigma_data, test_mode='flow_curve')

# Predict
sigma_pred = model.predict(gamma_dot)

# Get yield stress
tau_y = model.parameters.get_value("tau_y0")
print(f"Fitted yield stress: {tau_y:.1f} Pa")

Startup with Stress Overshoot

# Simulate startup
t = np.linspace(0, 50, 500)
gamma_dot = 1.0
t_wait = 100.0  # Waiting time before startup

strain, stress, fluidity = model.simulate_startup(t, gamma_dot, t_wait=t_wait)

# Analyze overshoot
sigma_max = np.max(stress)
sigma_ss = stress[-1]
overshoot_ratio = sigma_max / sigma_ss
print(f"Overshoot ratio: {overshoot_ratio:.2f}")

Creep Bifurcation

t = np.linspace(0, 1000, 500)

# Below yield - bounded strain
strain_below, _ = model.simulate_creep(t, sigma_applied=50.0)

# Above yield - unbounded flow
strain_above, _ = model.simulate_creep(t, sigma_applied=150.0)

# Plot shows bifurcation behavior

Normal Stress Predictions

The Fluidity-Saramito model tracks tensorial stress components (\(\tau_{xx}\), \(\tau_{yy}\), \(\tau_{xy}\)) internally via the upper-convected Maxwell constitutive equation, enabling first normal stress difference \(N_1 = \tau_{xx} - \tau_{yy}\) to emerge naturally during transient simulations (startup, LAOS). No public predict_normal_stresses() method is currently exposed; instead, access normal stress data through the tensorial stress output of simulate_startup() or simulate_laos().

LAOS Analysis

Large amplitude oscillatory shear (LAOS) provides a fingerprint of nonlinear material behavior through Lissajous curve shapes and harmonic decomposition.

Lissajous Curve Interpretation:

The Fluidity-Saramito model produces characteristic Lissajous shapes that reveal the interplay between elasticity, plasticity, and thixotropy:

Projection

Shape Feature

Physical Interpretation

Elastic (\(\sigma\) vs \(\gamma\))

Parallelogram deviation

Yielding, plastic flow above \(\tau_y\)

Elastic (\(\sigma\) vs \(\gamma\))

Enclosed area

Energy dissipation per cycle

Viscous (\(\sigma\) vs \(\dot{\gamma}\))

Bow-tie/self-intersection

Thixotropy, structure kinetics

Viscous (\(\sigma\) vs \(\dot{\gamma}\))

Asymmetric loops

Non-equilibrium fluidity dynamics

Elastic Projection ( \(\sigma\) vs \(\gamma\) ):

  • Linear viscoelastic: Ellipse with slope \(G'\) and area \(\propto G''\)

  • With yielding: Parallelogram-like corners where \(|\tau|\) crosses \(\tau_y\)

  • With thixotropy: Cycle-to-cycle evolution as fluidity equilibrates

  • Strain softening: Counterclockwise deviation (intracycle) indicates structure breakdown

Viscous Projection ( \(\sigma\) vs \(\dot{\gamma}\) ):

  • Linear viscoelastic: Ellipse with slope \(\eta'\) and area \(\propto \eta''\)

  • With thixotropy: Self-intersecting “bow-tie” or “figure-8” patterns

  • Interpretation: Self-intersection occurs when fluidity lags the deformation, causing stress to increase/decrease non-monotonically during each quarter cycle

  • Secondary loops: Indicate multiple structural relaxation timescales

Thixotropic Loop Signatures:

Pattern

Cause

Parameter Sensitivity

Self-intersecting viscous loops

\(f\) lags \(\dot{\gamma}\) by \(\pi/2\)

\(\tau_{\text{age}}\) vs \(1/\omega\) ratio

Cycle evolution (softening)

\(f\) not at steady state

Number of pre-cycles needed

Asymmetric peaks

Non-symmetric \(f(t)\) waveform

\(n_{\text{rej}} \neq 1\)

Harmonic distortion

Nonlinear plasticity

\(\gamma_0/\gamma_y\) ratio

Harmonic Extraction:

Fourier decomposition of the stress waveform provides quantitative nonlinearity measures:

from rheojax.models.fluidity import FluiditySaramitoLocal
import numpy as np

model = FluiditySaramitoLocal(coupling="minimal")
gamma_0 = 1.0  # 100% strain amplitude
omega = 1.0    # rad/s

# Simulate LAOS with sufficient cycles for steady state
t, strain, stress = model.simulate_laos(gamma_0, omega, n_cycles=10, n_points_per_cycle=100)

# Extract harmonics from last cycle (steady state)
last_cycle_idx = -100  # Last 100 points = 1 cycle
stress_cycle = stress[last_cycle_idx:]

# FFT for harmonic amplitudes
fft_stress = np.fft.fft(stress_cycle)
harmonics = np.abs(fft_stress[:5]) / len(stress_cycle) * 2

print(f"1st harmonic (G*): {harmonics[1]:.1f} Pa")
print(f"3rd harmonic (I_3): {harmonics[3]:.1f} Pa")
print(f"I_3/I_1 ratio: {harmonics[3]/harmonics[1]:.4f}")

# Chebyshev coefficients for LAOS analysis
# e_1 (strain-stiffening) and v_1 (shear-thickening)
from numpy.polynomial import chebyshev
gamma_norm = strain[last_cycle_idx:] / gamma_0
cheb_coeffs = chebyshev.chebfit(gamma_norm, stress_cycle, deg=5)

Standard LAOS Example:

from rheojax.models.fluidity import FluiditySaramitoLocal
import numpy as np
import matplotlib.pyplot as plt

# Large amplitude oscillatory shear
model = FluiditySaramitoLocal(coupling="minimal")

gamma_0 = 1.0  # 100% strain
omega = 1.0    # rad/s

# Simulate LAOS
t, strain, stress = model.simulate_laos(gamma_0, omega, n_cycles=3)

# Plot Lissajous curves
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Elastic projection (σ vs γ)
axes[0].plot(strain, stress)
axes[0].set_xlabel('Strain γ')
axes[0].set_ylabel('Stress σ (Pa)')
axes[0].set_title('Elastic Lissajous (σ vs γ)')
axes[0].grid(True)

# Viscous projection (σ vs γ̇)
gamma_dot = gamma_0 * omega * np.cos(omega * t)
axes[1].plot(gamma_dot, stress)
axes[1].set_xlabel('Strain rate γ̇ (1/s)')
axes[1].set_ylabel('Stress σ (Pa)')
axes[1].set_title('Viscous Lissajous (σ vs γ̇)')
axes[1].grid(True)

plt.tight_layout()

Nonlocal Model with Shear Banding

from rheojax.models.fluidity import FluiditySaramitoNonlocal

# Create nonlocal model
model = FluiditySaramitoNonlocal(
    coupling="minimal",
    N_y=51,      # Grid points
    H=1e-3,      # Gap width (m)
    xi=1e-5,     # Cooperativity length (m)
)

# Simulate startup
t = np.linspace(0, 50, 200)
_, sigma, f_field = model.simulate_startup(t, gamma_dot=0.1)

# Check for shear banding
is_banded, cv, ratio = model.detect_shear_bands()
print(f"Shear banding: {is_banded}, CV={cv:.2f}, ratio={ratio:.1f}")

# Get detailed metrics
metrics = model.get_banding_metrics()
print(f"Band fraction: {metrics['band_fraction']:.2f}")

Bayesian Inference

# Fit with NLSQ first (warm-start)
model.fit(gamma_dot, sigma, test_mode='flow_curve')

# Bayesian inference
result = model.fit_bayesian(
    gamma_dot, sigma,
    test_mode='flow_curve',
    num_warmup=1000,
    num_samples=2000,
    num_chains=4,  # Production-ready diagnostics
    seed=42,       # Reproducibility
)

# Check diagnostics
# R-hat should be < 1.01, ESS > 400

# Plot with ArviZ
from rheojax.pipeline.bayesian import BayesianPipeline

pipeline = BayesianPipeline()
pipeline._idata = result.idata
pipeline.plot_trace()
pipeline.plot_pair(divergences=True)

Comparison with Existing Models

vs FluidityLocal

Feature

FluidityLocal

FluiditySaramitoLocal

Stress tensor

Scalar \(\sigma\)

Tensorial [\(\tau_{xx}\), \(\tau_{yy}\), \(\tau_{xy}\)]

Normal stresses

No

Yes (\(N_1\) from UCM)

Yield criterion

None (implicit)

Von Mises (explicit)

Elastic effects

Maxwell-like

Upper-convected Maxwell

Parameters

9

10-12

The Saramito model is preferred when:

  • Normal stresses (\(N_1\)) are important

  • Tensorial stress state is needed

  • True yield criterion is required

vs Standard Saramito (no fluidity)

The fluidity extension adds:

  • Thixotropic time dependence

  • Aging/rejuvenation dynamics

  • Shear banding capability (nonlocal)

Without fluidity, the standard Saramito model has constant relaxation time and cannot capture thixotropic behaviors.

vs Other Yield-Stress and Thixotropic Models

The Fluidity-Saramito model occupies a unique position in the landscape of yield-stress and thixotropic constitutive models. This comprehensive comparison helps practitioners select the appropriate model.

Feature Comparison Table:

Feature

Fluidity-Saramito

STZ

HL Trap

Giesekus

DMT

Yield stress

Explicit (\(\tau_y\))

Emergent

Emergent

None

Explicit

Aging mechanism

Explicit (\(f\) kinetics)

Via \(\chi\) (effective temp)

Via \(n(E)\) trap distribution

None

Via \(\lambda\) kinetics

Stress tensor

Full tensorial

Scalar (typically)

Scalar

Full tensorial

Scalar

Shear banding

With \(D_f\) (nonlocal)

Needs spatial extension

Needs spatial extension

No

With nonlocal extension

Normal stresses

Yes (\(N_1\))

Limited

No

Yes (\(N_1, N_2\))

No

# Parameters

10-12

~7

~7

~4

~8

Physical basis

Phenomenological + UCM

Statistical mechanics

Mean-field kinetic

Molecular network

Structural kinetics

Model Selection Guide:

When to use Fluidity-Saramito:

  • Tensorial stress required: Complex flow geometries, extensional components, normal stress predictions for rod climbing, die swell

  • Explicit yield criterion: Need to track Von Mises stress vs \(\tau_y\)

  • Clear separation of effects: Distinct elastic, plastic, and structural contributions identifiable from parameters

  • Dense colloids with both aging and elasticity: Emulsions, polymer gels, cosmetics where both \(G\) and thixotropy matter

  • Validation against velocimetry: Nonlocal variant directly predicts \(f(y)\) profiles comparable to PIV/USV data

When to use alternatives:

Model

Use When

STZ

Metallic glasses, granular media; need statistical mechanics foundation; yield emerges from disorder physics

HL Trap

Colloidal glasses near jamming; energy landscape interpretation; distribution of relaxation times

Giesekus

Polymer melts/solutions without yield; need \(N_2\) prediction; fewer parameters, no thixotropy

DMT

Similar physics to Fluidity-Saramito but with scalar stress and structure parameter \(\lambda \in [0,1]\); simpler tensorial extension

Complexity vs. Capability Trade-off:

Capability (features modeled)
      ▲
      │                                    ┌─────────────────────┐
      │                                    │ Fluidity-Saramito   │
      │                              ┌─────┤ Nonlocal            │
      │                              │     └─────────────────────┘
      │                    ┌─────────┴───┐
      │                    │ Fluidity-   │
      │                    │ Saramito    │
      │          ┌─────────┤ Local       │
      │          │         └─────────────┘
      │    ┌─────┴────┐     ┌──────┐
      │    │ DMT/     │     │ STZ  │
      │    │ Fluidity │     └──────┘
      │    │ Local    │  ┌───────┐
      │    └──────────┘  │ HL    │
      │                  └───────┘
      │  ┌──────────┐
      │  │ HB/      │
      │  │ Bingham  │
      │  └──────────┘
      └───────────────────────────────────────────────────────────► Complexity
                                                            (# parameters)

References

API Reference

class rheojax.models.fluidity.saramito.FluiditySaramitoLocal(coupling='minimal')[source]

Local (0D) Fluidity-Saramito Model for elastoviscoplastic fluids.

Implements a homogeneous Saramito EVP model where the material state is characterized by tensorial stress [τ_xx, τ_yy, τ_xy] and scalar fluidity f(t) that evolve via coupled ODEs.

The constitutive equation (Saramito 2009): λ(∇̂τ) + α(τ)τ = 2η_p D

where:

  • λ = 1/f: Fluidity-dependent relaxation time

  • ∇̂τ: Upper-convected derivative

  • α = max(0, 1 - τ_y/|τ|): Von Mises plasticity

  • η_p = G/f: Polymeric viscosity

Fluidity evolves via: df/dt = (f_age - f)/t_a + b|γ̇|^n(f_flow - f)

Parameters:

coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution: - “minimal”: τ_y = τ_y0 (constant) - “full”: τ_y(f) = τ_y0 + a_y/f^m (aging yield stress)

Examples

Basic flow curve fitting:

>>> model = FluiditySaramitoLocal(coupling="minimal")
>>> model.fit(gamma_dot, sigma, test_mode="flow_curve")
>>> sigma_pred = model.predict(gamma_dot)

Startup with stress overshoot:

>>> model = FluiditySaramitoLocal()
>>> model.fit(t, sigma, test_mode="startup", gamma_dot=1.0)
>>> strain, stress, f = model.simulate_startup(t, gamma_dot=1.0)

Bayesian inference:

>>> result = model.fit_bayesian(gamma_dot, sigma, test_mode="flow_curve",
...                             num_warmup=1000, num_samples=2000)
>>> intervals = model.get_credible_intervals(result.posterior_samples)
__init__(coupling='minimal')[source]

Initialize Local Fluidity-Saramito Model.

Parameters:

coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution

simulate_startup(t, gamma_dot, t_wait=0.0)[source]

Simulate startup response with full trajectory.

Parameters:

  • t (ndarray) – Time array (s)

  • gamma_dot (float) – Applied shear rate (1/s)

  • t_wait (float) – Waiting time before startup (s)

Return type:

tuple[ndarray, ndarray, ndarray]

Returns:

  • strain (np.ndarray) – Accumulated strain γ(t)

  • stress (np.ndarray) – Shear stress τ_xy(t) (Pa)

  • fluidity (np.ndarray) – Fluidity f(t) (1/(Pa·s))

simulate_creep(t, sigma_applied, t_wait=0.0)[source]

Simulate creep response.

Parameters:

  • t (ndarray) – Time array (s)

  • sigma_applied (float) – Applied stress (Pa)

  • t_wait (float) – Waiting time before creep (s)

Return type:

tuple[ndarray, ndarray]

Returns:

  • strain (np.ndarray) – Accumulated strain γ(t)

  • fluidity (np.ndarray) – Fluidity f(t) (1/(Pa·s))

simulate_laos(gamma_0, omega, n_cycles=3, n_points_per_cycle=256)[source]

Simulate LAOS response.

Parameters:

  • gamma_0 (float) – Strain amplitude

  • omega (float) – Angular frequency (rad/s)

  • n_cycles (int) – Number of oscillation cycles

  • n_points_per_cycle (int) – Points per cycle

Return type:

tuple[ndarray, ndarray, ndarray]

Returns:

  • t (np.ndarray) – Time array (s)

  • strain (np.ndarray) – Strain array

  • stress (np.ndarray) – Stress array (Pa)

extract_harmonics(stress, n_points_per_cycle=256)[source]

Extract Fourier harmonics from LAOS stress response.

Parameters:

  • stress (ndarray) – Stress array from simulate_laos

  • n_points_per_cycle (int) – Points per cycle

Returns:

Dictionary with I_1, I_3, I_5 amplitudes and ratios

Return type:

dict

model_function(X, params, test_mode=None, **kwargs)[source]

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode. Accepts protocol-specific kwargs (gamma_dot, sigma_applied, etc.).

Parameters:

  • X (array-like) – Independent variable

  • params (array-like) – Parameter values in order

  • test_mode (str, optional) – Override stored test mode

  • **kwargs – Protocol-specific arguments (gamma_dot, sigma_applied, gamma_0, omega)

Returns:

Predicted response

Return type:

jnp.ndarray

get_overshoot_ratio(gamma_dot, t_max=100.0, n_points=1000)[source]

Compute stress overshoot ratio σ_max / σ_ss.

Parameters:

  • gamma_dot (float) – Shear rate (1/s)

  • t_max (float) – Maximum simulation time (s)

  • n_points (int) – Number of time points

Returns:

Overshoot ratio (1.0 = no overshoot)

Return type:

float

get_critical_stress()[source]

Get critical stress for creep bifurcation.

Returns:

Critical stress σ_c ≈ τ_y (Pa)

Return type:

float

__repr__()

Return string representation.

Return type:

str

property deborah_number: float | None

Get Deborah number De = λ * ω (for oscillatory tests).

Returns:

Deborah number if omega is set, else None

Return type:

float or None

fit(X, y, method='nlsq', check_compatibility=False, use_log_residuals=None, use_multi_start=None, n_starts=5, perturb_factor=0.3, deformation_mode=None, poisson_ratio=0.5, auto_init=False, return_result=False, check_physics=False, uncertainty=None, **kwargs)

Fit the model to data using NLSQ optimization.

This method uses NLSQ (GPU-accelerated nonlinear least squares) by default for fast point estimation. The optimization result is stored for potential warm-starting of Bayesian inference.

For very wide frequency ranges (>10 decades), multi-start optimization is automatically enabled to escape local minima.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • method (str) – Optimization method (‘nlsq’ by default for compatibility)

  • check_compatibility (bool) – Whether to check model-data compatibility before fitting. If True, warns when model may not be appropriate for data. Default is False for backward compatibility.

  • use_log_residuals (bool | None) – Whether to use log-space residuals for fitting. Recommended for wide frequency ranges (>8 decades) to prevent optimizer bias. If None (default), automatically detected based on data range. Explicit True/False overrides auto-detection.

  • use_multi_start (bool | None) – Whether to use multi-start optimization to escape local minima. Recommended for very wide ranges (>10 decades). If None (default), automatically enabled for >10 decades.

  • n_starts (int) – Number of random starts for multi-start optimization (default: 5)

  • perturb_factor (float) – Perturbation magnitude for multi-start random starts (default: 0.3). Parameters are perturbed by ± perturb_factor * (value or range). Larger values (0.7-0.9) explore wider parameter space.

  • auto_init (bool) – If True, calls auto_p0() to estimate initial parameters from data before running the optimizer (default: False).

  • return_result (bool) – If True, returns a FitResult instead of self. This intentionally breaks method chaining for workflows that need structured result objects (default: False).

  • check_physics (bool) – If True, runs post-fit physics validation and emits RheoJaxPhysicsWarning for any violations (default: False).

  • uncertainty (str | None) – Post-fit uncertainty method. "hessian" for fast Cramér-Rao bounds, "bootstrap" for residual bootstrap CIs, or None to skip (default: None).

  • **kwargs – Additional fitting options passed to _fit()

Return type:

BaseModel | Any

Returns:

self for method chaining (default), or FitResult if return_result=True.

Example

>>> model = Maxwell()
>>> model.fit(t, G_data)  # Uses NLSQ by default
>>> model.fit(t, G_data, method='nlsq', max_iter=1000)
>>> model.fit(t, G_data, check_compatibility=True)  # Check compatibility
>>> model.fit(omega, G_star, use_log_residuals=True)  # Force log-residuals
>>> model.fit(mastercurve, None, use_multi_start=True, n_starts=10)  # Multi-start
>>> result = model.fit(t, G_data, return_result=True)  # Structured result
>>> result = model.fit(t, G_data, auto_init=True, check_physics=True,
...                    return_result=True)  # Full pipeline
fit_bayesian(X, y=None, num_warmup=1000, num_samples=2000, num_chains=4, initial_values=None, test_mode=None, seed=None, deformation_mode=None, poisson_ratio=0.5, **nuts_kwargs)

Perform Bayesian inference using NumPyro NUTS sampler.

This method delegates to BayesianMixin.fit_bayesian() to run NUTS sampling for Bayesian parameter estimation. If initial_values is not provided and the model has been previously fitted with fit(), the NLSQ point estimates are automatically used for warm-starting.

Multi-chain sampling is enabled by default (num_chains=4) to provide reliable convergence diagnostics (R-hat, ESS) and parallel execution on multi-GPU systems.

Parameters:

  • X (numpy.typing.ArrayLike) – Independent variable data (input features) or RheoData object

  • y (Optional[numpy.typing.ArrayLike]) – Dependent variable data (observations to fit). If X is RheoData, y is ignored and extracted from X.

  • num_warmup (int) – Number of warmup/burn-in iterations (default: 1000)

  • num_samples (int) – Number of posterior samples per chain (default: 2000)

  • num_chains (int) – Number of MCMC chains (default: 4). Multiple chains enable proper R-hat computation and parallel execution. Chain method is auto-selected: ‘parallel’ on multi-GPU, ‘vectorized’ on single GPU/CPU.

  • initial_values (dict[str, float] | None) – Optional dict of initial parameter values for warm-start. If None and model is fitted, uses NLSQ estimates.

  • test_mode (str | None) – Explicit test mode (e.g., ‘relaxation’, ‘creep’, ‘oscillation’). If None, inferred from RheoData.metadata[‘test_mode’] or defaults to ‘relaxation’. Overrides RheoData metadata if provided.

  • seed (int | None) – Random seed for reproducibility. If None, uses seed=0 for deterministic results. Set to different values for independent runs.

  • **nuts_kwargs – Additional arguments passed to NUTS sampler (e.g., target_accept_prob, chain_method)

Return type:

BayesianResult

Returns:

BayesianResult containing posterior samples, summary statistics, and convergence diagnostics (R-hat, ESS, divergences)

Example

>>> model = Maxwell()
>>> # Warm-start from NLSQ with explicit mode
>>> model.fit(t, G_data, test_mode='relaxation')  # NLSQ optimization
>>> result = model.fit_bayesian(t, G_data, test_mode='relaxation')
>>>
>>> # RheoData with embedded mode (recommended)
>>> rheo_data = RheoData(x=omega, y=G_star, metadata={'test_mode': 'oscillation'})
>>> result = model.fit_bayesian(rheo_data)
>>>
>>> # Or provide explicit initial values
>>> result = model.fit_bayesian(
...     t, G_data,
...     initial_values={'G0': 1e5, 'eta': 1e3},
...     test_mode='creep'
... )
fit_predict(X, y, **kwargs)

Fit model and return predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • **kwargs – Additional fitting options

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions on training data

classmethod from_dict(data)

Create model from dictionary.

Parameters:

data (dict[str, Any]) – Dictionary representation

Return type:

BaseModel

Returns:

Model instance

get_bayesian_result()

Get stored Bayesian inference result.

Return type:

BayesianResult | None

Returns:

BayesianResult from fit_bayesian(), or None if not run

Example

>>> model.fit_bayesian(t, G_data)
>>> result = model.get_bayesian_result()
>>> print(result.diagnostics['r_hat'])
get_credible_intervals(posterior_samples, credibility=0.95)

Compute highest density intervals (HDI) for posterior samples.

Computes the highest posterior density interval for each parameter, which is the shortest interval containing the specified probability mass. This is preferred over equal-tailed intervals for skewed distributions.

Parameters:

  • posterior_samples (dict[str, ndarray]) – Dictionary mapping parameter names to posterior sample arrays (from BayesianResult.posterior_samples)

  • credibility (float) – Probability mass to include in interval (default: 0.95) Common values: 0.68 (1 sigma), 0.95 (2 sigma), 0.997 (3 sigma)

Return type:

dict[str, tuple[float, float]]

Returns:

Dictionary mapping parameter names to (lower, upper) credible interval tuples. All values are float64.

Example

>>> result = model.fit_bayesian(X, y)
>>> intervals_95 = model.get_credible_intervals(
...     result.posterior_samples, credibility=0.95
... )
>>> print(f"95% CI for a: {intervals_95['a']}")
get_effective_yield_stress(f)

Get effective yield stress based on coupling mode.

Parameters:

f (float) – Current fluidity (1/(Pa·s))

Returns:

Effective yield stress τ_y(f) (Pa)

Return type:

float

get_initial_fluidity(t_wait=0.0)

Get initial fluidity value based on waiting time.

For aged samples (t_wait >> t_a): f → f_age For fresh samples (t_wait = 0): f can start at f_flow

Parameters:

t_wait (float) – Waiting time before measurement (s)

Returns:

Initial fluidity f_0 (1/(Pa·s))

Return type:

float

get_nlsq_result()

Get stored NLSQ optimization result.

Returns:

OptimizationResult from NLSQ fit, or None if not fitted

Example

>>> model.fit(t, G_data)
>>> result = model.get_nlsq_result()
>>> if result:
...     print(f"Converged: {result.success}")
get_parameter_dict()

Get all parameters as a dictionary.

Returns:

Dictionary of parameter name → value

Return type:

dict

get_parameter_uncertainties()

Get standard errors for fitted parameters from NLSQ covariance.

Returns:

std_error}, or None if covariance unavailable

Return type:

dict of {param_name

get_params(deep=True)

Get model parameters.

Parameters:

deep (bool) – If True, return parameters of sub-objects

Return type:

dict[str, Any]

Returns:

Dictionary of parameter names and values

initialize_from_flow_curve(gamma_dot, sigma)

Initialize parameters from flow curve data.

Smart initialization strategy: 1. τ_y0 from low-rate plateau (stress intercept) 2. K_HB, n_HB from high-rate power-law region 3. G estimated from stress overshoot if available 4. Fluidity parameters from viscosity scaling

Parameters:

  • gamma_dot (ndarray) – Shear rate data (1/s)

  • sigma (ndarray) – Shear stress data (Pa)

Return type:

None

property pcov_

Parameter covariance matrix from NLSQ fit.

Returns:

ndarray of shape (n_params, n_params), or None if not fitted

property popt_

Optimal parameter values from NLSQ fit.

Returns:

ndarray of shape (n_params,), or None if not fitted

precompile(test_mode='relaxation', X=None, y=None)

Precompile NLSQ residual functions to eliminate JIT cold-start.

Triggers JIT compilation by running a minimal fit (max_iter=1) with dummy data. The model parameters are reset afterwards so the model is left in its original state.

This is useful for interactive sessions or benchmarks where the ~870ms first-fit JIT overhead should be excluded.

Parameters:

  • test_mode (str) – Test mode to precompile for (default: ‘relaxation’).

  • X (Optional[numpy.typing.ArrayLike]) – Optional input data for shape inference. If None, uses a 10-point logspace array.

  • y (Optional[numpy.typing.ArrayLike]) – Optional output data. If None, generates ones matching X.

Return type:

float

Returns:

Compilation time in seconds.

Example

>>> model = Maxwell()
>>> t = model.precompile(test_mode='relaxation')
>>> print(f"Compiled in {t:.2f}s")
>>> model.fit(X, y)  # No JIT overhead
precompile_bayesian(X=None, y=None, test_mode=None, num_chains=4)

Precompile NUTS kernel to eliminate JIT overhead in subsequent calls.

Triggers JIT compilation of the NumPyro model by running a minimal sampling (1 warmup, 1 sample). This caches the compiled kernel so that subsequent fit_bayesian() calls are 2-5x faster.

Parameters:

  • X (ndarray | RheoData | None) – Sample input data for determining array shapes. If None, uses a default 10-point linspace [0.01, 100].

  • y (ndarray | None) – Sample output data. If None, generates dummy data.

  • test_mode (str | TestModeEnum | None) – Test mode to precompile for. If None, defaults to ‘relaxation’.

Returns:

Compilation time in seconds.

Return type:

float

Example

>>> model = Maxwell()
>>> compile_time = model.precompile_bayesian(test_mode='relaxation')
>>> print(f"Compiled in {compile_time:.1f}s")
>>> # Now fit_bayesian() will be faster
>>> result = model.fit_bayesian(X, y)  # No compilation overhead
predict(X, test_mode=None, deformation_mode=None, poisson_ratio=None, **kwargs)

Make predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • test_mode (str | None) – Optional test mode (‘oscillation’, ‘relaxation’, ‘creep’, ‘flow’). If provided, sets model’s test_mode before prediction. Useful for data generation without fitting.

  • deformation_mode (str | DeformationMode | None) – Optional deformation mode for output conversion. If None, uses the mode stored from fit(). If tensile, converts G* predictions to E* space.

  • poisson_ratio (float | None) – Poisson’s ratio for conversion. If None, uses value stored from fit() (default 0.5).

  • **kwargs – Additional arguments passed to the internal _predict method.

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions (in E* space if deformation_mode is tensile)

predict_normal_stresses(gamma_dot)

Predict first and second normal stress differences.

The Saramito model predicts non-zero N₁ from its upper-convected Maxwell foundation.

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)

Return type:

tuple[ndarray, ndarray]

Returns:

  • N1 (np.ndarray) – First normal stress difference τ_xx - τ_yy (Pa)

  • N2 (np.ndarray) – Second normal stress difference τ_yy - τ_zz (Pa) Note: N2 = 0 for simple UCM models

property relaxation_time: float

Get characteristic relaxation time λ = 1/(G*f_age).

This is the relaxation time in the aged (rest) state.

Returns:

Relaxation time (s)

Return type:

float

sample_prior(num_samples=1000, seed=None)

Sample from prior distributions over parameter bounds.

Samples from uniform prior distributions defined by parameter bounds. This is useful for prior predictive checks and understanding the prior’s influence on the posterior.

Parameters:

  • num_samples (int) – Number of samples to draw from prior (default: 1000)

  • seed (int | None) – Random seed for reproducibility (default: None)

Return type:

dict[str, ndarray]

Returns:

Dictionary mapping parameter names to arrays of prior samples. Each array has shape (num_samples,) and dtype float64.

Raises:

Example

>>> model = MyModel()
>>> prior_samples = model.sample_prior(num_samples=500, seed=42)
>>> print(prior_samples["a"].shape)  # (500,)
score(X, y)

Compute model score (R² by default).

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – True target values

Return type:

float

Returns:

Model score (R² coefficient)

set_params(**params)

Set model parameters.

Parameters:

**params – Parameter names and values

Return type:

BaseModel

Returns:

self for method chaining

to_dict()

Serialize model to dictionary.

Return type:

dict[str, Any]

Returns:

Dictionary representation of model

property weissenberg_number: float | None

Get Weissenberg number Wi = λ * γ̇ (for steady/startup tests).

Returns:

Weissenberg number if gamma_dot is set, else None

Return type:

float or None

parameters: ParameterSet
class rheojax.models.fluidity.saramito.FluiditySaramitoNonlocal(coupling='minimal', N_y=51, H=0.001, xi=1e-05)[source]

Nonlocal (1D) Fluidity-Saramito Model with spatial diffusion.

Implements a spatially-resolved Saramito EVP model where fluidity varies across a Couette gap and can form shear bands.

The fluidity evolution includes a diffusion term: ∂f/∂t = (f_loc - f)/t_a + D_f * ∇²f

where D_f = ξ²/t_a is the fluidity diffusivity and ξ is the cooperativity length (interface width parameter).

Parameters:

  • coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution

  • N_y (int) – Number of spatial grid points

  • H (float) – Gap width (m)

  • xi (float) – Cooperativity length (m)

Notes

The model solves a coupled PDE system for [Σ, f(y)] where Σ is the bulk (gap-averaged) stress. In Couette geometry, stress is approximately uniform across the gap, enabling this simplification.

Shear bands appear when the fluidity profile develops large gradients, with high-fluidity (flowing) bands coexisting with low-fluidity (jammed) regions.

Examples

Basic flow curve with banding check:

>>> model = FluiditySaramitoNonlocal(N_y=51, H=1e-3, xi=1e-5)
>>> model.fit(gamma_dot, sigma, test_mode="flow_curve")
>>> is_banded, cv, ratio = model.detect_shear_bands()
>>> print(f"Shear banding detected: {is_banded}")

Startup transient with spatial profile:

>>> model = FluiditySaramitoNonlocal()
>>> t, sigma, f_field = model.simulate_startup(t, gamma_dot=1.0)
>>> model.plot_fluidity_profile()  # Shows spatial variation
__init__(coupling='minimal', N_y=51, H=0.001, xi=1e-05)[source]

Initialize Nonlocal Fluidity-Saramito Model.

Parameters:

  • coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution

  • N_y (int) – Number of spatial grid points

  • H (float) – Gap width (m)

  • xi (float) – Cooperativity length (m)

simulate_startup(t, gamma_dot)[source]

Simulate startup with spatial resolution.

Parameters:

  • t (ndarray) – Time array (s)

  • gamma_dot (float) – Applied shear rate (1/s)

Return type:

tuple[ndarray, ndarray, ndarray]

Returns:

  • t (np.ndarray) – Time array

  • sigma (np.ndarray) – Bulk stress (Pa)

  • f_field (np.ndarray) – Final fluidity profile, shape (N_y,)

simulate_creep(t, sigma_applied)[source]

Simulate creep with spatial resolution.

Parameters:

  • t (ndarray) – Time array (s)

  • sigma_applied (float) – Applied stress (Pa)

Return type:

tuple[ndarray, ndarray]

Returns:

  • gamma (np.ndarray) – Bulk strain

  • f_field (np.ndarray) – Final fluidity profile

detect_shear_bands(f_profile=None, threshold=0.3)[source]

Detect shear banding from fluidity profile.

Parameters:

  • f_profile (ndarray | None) – Fluidity field. If None, uses stored profile.

  • threshold (float) – CV threshold for banding detection

Return type:

tuple[bool, float, float]

Returns:

  • is_banded (bool) – True if shear banding detected

  • cv (float) – Coefficient of variation

  • ratio (float) – Max/min fluidity ratio

get_banding_metrics(f_profile=None)[source]

Get detailed shear banding metrics.

Parameters:

f_profile (ndarray | None) – Fluidity field. If None, uses stored profile.

Returns:

Metrics including cv, ratio, band_width, etc.

Return type:

dict[str, float]

model_function(X, params, test_mode=None, **kwargs)[source]

NumPyro/BayesianMixin model function.

Accepts protocol-specific kwargs (gamma_dot, sigma_applied, etc.).

Parameters:

  • X (array-like) – Independent variable

  • params (array-like) – Parameter values

  • test_mode (str, optional) – Override stored test mode

  • **kwargs – Protocol-specific arguments (gamma_dot, sigma_applied)

Returns:

Predicted response

Return type:

jnp.ndarray

property y_grid: ndarray

Get spatial grid across gap.

Returns:

Position array (m)

Return type:

np.ndarray

__repr__()[source]

Return string representation.

Return type:

str

property deborah_number: float | None

Get Deborah number De = λ * ω (for oscillatory tests).

Returns:

Deborah number if omega is set, else None

Return type:

float or None

fit(X, y, method='nlsq', check_compatibility=False, use_log_residuals=None, use_multi_start=None, n_starts=5, perturb_factor=0.3, deformation_mode=None, poisson_ratio=0.5, auto_init=False, return_result=False, check_physics=False, uncertainty=None, **kwargs)

Fit the model to data using NLSQ optimization.

This method uses NLSQ (GPU-accelerated nonlinear least squares) by default for fast point estimation. The optimization result is stored for potential warm-starting of Bayesian inference.

For very wide frequency ranges (>10 decades), multi-start optimization is automatically enabled to escape local minima.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • method (str) – Optimization method (‘nlsq’ by default for compatibility)

  • check_compatibility (bool) – Whether to check model-data compatibility before fitting. If True, warns when model may not be appropriate for data. Default is False for backward compatibility.

  • use_log_residuals (bool | None) – Whether to use log-space residuals for fitting. Recommended for wide frequency ranges (>8 decades) to prevent optimizer bias. If None (default), automatically detected based on data range. Explicit True/False overrides auto-detection.

  • use_multi_start (bool | None) – Whether to use multi-start optimization to escape local minima. Recommended for very wide ranges (>10 decades). If None (default), automatically enabled for >10 decades.

  • n_starts (int) – Number of random starts for multi-start optimization (default: 5)

  • perturb_factor (float) – Perturbation magnitude for multi-start random starts (default: 0.3). Parameters are perturbed by ± perturb_factor * (value or range). Larger values (0.7-0.9) explore wider parameter space.

  • auto_init (bool) – If True, calls auto_p0() to estimate initial parameters from data before running the optimizer (default: False).

  • return_result (bool) – If True, returns a FitResult instead of self. This intentionally breaks method chaining for workflows that need structured result objects (default: False).

  • check_physics (bool) – If True, runs post-fit physics validation and emits RheoJaxPhysicsWarning for any violations (default: False).

  • uncertainty (str | None) – Post-fit uncertainty method. "hessian" for fast Cramér-Rao bounds, "bootstrap" for residual bootstrap CIs, or None to skip (default: None).

  • **kwargs – Additional fitting options passed to _fit()

Return type:

BaseModel | Any

Returns:

self for method chaining (default), or FitResult if return_result=True.

Example

>>> model = Maxwell()
>>> model.fit(t, G_data)  # Uses NLSQ by default
>>> model.fit(t, G_data, method='nlsq', max_iter=1000)
>>> model.fit(t, G_data, check_compatibility=True)  # Check compatibility
>>> model.fit(omega, G_star, use_log_residuals=True)  # Force log-residuals
>>> model.fit(mastercurve, None, use_multi_start=True, n_starts=10)  # Multi-start
>>> result = model.fit(t, G_data, return_result=True)  # Structured result
>>> result = model.fit(t, G_data, auto_init=True, check_physics=True,
...                    return_result=True)  # Full pipeline
fit_bayesian(X, y=None, num_warmup=1000, num_samples=2000, num_chains=4, initial_values=None, test_mode=None, seed=None, deformation_mode=None, poisson_ratio=0.5, **nuts_kwargs)

Perform Bayesian inference using NumPyro NUTS sampler.

This method delegates to BayesianMixin.fit_bayesian() to run NUTS sampling for Bayesian parameter estimation. If initial_values is not provided and the model has been previously fitted with fit(), the NLSQ point estimates are automatically used for warm-starting.

Multi-chain sampling is enabled by default (num_chains=4) to provide reliable convergence diagnostics (R-hat, ESS) and parallel execution on multi-GPU systems.

Parameters:

  • X (numpy.typing.ArrayLike) – Independent variable data (input features) or RheoData object

  • y (Optional[numpy.typing.ArrayLike]) – Dependent variable data (observations to fit). If X is RheoData, y is ignored and extracted from X.

  • num_warmup (int) – Number of warmup/burn-in iterations (default: 1000)

  • num_samples (int) – Number of posterior samples per chain (default: 2000)

  • num_chains (int) – Number of MCMC chains (default: 4). Multiple chains enable proper R-hat computation and parallel execution. Chain method is auto-selected: ‘parallel’ on multi-GPU, ‘vectorized’ on single GPU/CPU.

  • initial_values (dict[str, float] | None) – Optional dict of initial parameter values for warm-start. If None and model is fitted, uses NLSQ estimates.

  • test_mode (str | None) – Explicit test mode (e.g., ‘relaxation’, ‘creep’, ‘oscillation’). If None, inferred from RheoData.metadata[‘test_mode’] or defaults to ‘relaxation’. Overrides RheoData metadata if provided.

  • seed (int | None) – Random seed for reproducibility. If None, uses seed=0 for deterministic results. Set to different values for independent runs.

  • **nuts_kwargs – Additional arguments passed to NUTS sampler (e.g., target_accept_prob, chain_method)

Return type:

BayesianResult

Returns:

BayesianResult containing posterior samples, summary statistics, and convergence diagnostics (R-hat, ESS, divergences)

Example

>>> model = Maxwell()
>>> # Warm-start from NLSQ with explicit mode
>>> model.fit(t, G_data, test_mode='relaxation')  # NLSQ optimization
>>> result = model.fit_bayesian(t, G_data, test_mode='relaxation')
>>>
>>> # RheoData with embedded mode (recommended)
>>> rheo_data = RheoData(x=omega, y=G_star, metadata={'test_mode': 'oscillation'})
>>> result = model.fit_bayesian(rheo_data)
>>>
>>> # Or provide explicit initial values
>>> result = model.fit_bayesian(
...     t, G_data,
...     initial_values={'G0': 1e5, 'eta': 1e3},
...     test_mode='creep'
... )
fit_predict(X, y, **kwargs)

Fit model and return predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • **kwargs – Additional fitting options

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions on training data

classmethod from_dict(data)

Create model from dictionary.

Parameters:

data (dict[str, Any]) – Dictionary representation

Return type:

BaseModel

Returns:

Model instance

get_bayesian_result()

Get stored Bayesian inference result.

Return type:

BayesianResult | None

Returns:

BayesianResult from fit_bayesian(), or None if not run

Example

>>> model.fit_bayesian(t, G_data)
>>> result = model.get_bayesian_result()
>>> print(result.diagnostics['r_hat'])
get_credible_intervals(posterior_samples, credibility=0.95)

Compute highest density intervals (HDI) for posterior samples.

Computes the highest posterior density interval for each parameter, which is the shortest interval containing the specified probability mass. This is preferred over equal-tailed intervals for skewed distributions.

Parameters:

  • posterior_samples (dict[str, ndarray]) – Dictionary mapping parameter names to posterior sample arrays (from BayesianResult.posterior_samples)

  • credibility (float) – Probability mass to include in interval (default: 0.95) Common values: 0.68 (1 sigma), 0.95 (2 sigma), 0.997 (3 sigma)

Return type:

dict[str, tuple[float, float]]

Returns:

Dictionary mapping parameter names to (lower, upper) credible interval tuples. All values are float64.

Example

>>> result = model.fit_bayesian(X, y)
>>> intervals_95 = model.get_credible_intervals(
...     result.posterior_samples, credibility=0.95
... )
>>> print(f"95% CI for a: {intervals_95['a']}")
get_effective_yield_stress(f)

Get effective yield stress based on coupling mode.

Parameters:

f (float) – Current fluidity (1/(Pa·s))

Returns:

Effective yield stress τ_y(f) (Pa)

Return type:

float

get_initial_fluidity(t_wait=0.0)

Get initial fluidity value based on waiting time.

For aged samples (t_wait >> t_a): f → f_age For fresh samples (t_wait = 0): f can start at f_flow

Parameters:

t_wait (float) – Waiting time before measurement (s)

Returns:

Initial fluidity f_0 (1/(Pa·s))

Return type:

float

get_nlsq_result()

Get stored NLSQ optimization result.

Returns:

OptimizationResult from NLSQ fit, or None if not fitted

Example

>>> model.fit(t, G_data)
>>> result = model.get_nlsq_result()
>>> if result:
...     print(f"Converged: {result.success}")
get_parameter_dict()

Get all parameters as a dictionary.

Returns:

Dictionary of parameter name → value

Return type:

dict

get_parameter_uncertainties()

Get standard errors for fitted parameters from NLSQ covariance.

Returns:

std_error}, or None if covariance unavailable

Return type:

dict of {param_name

get_params(deep=True)

Get model parameters.

Parameters:

deep (bool) – If True, return parameters of sub-objects

Return type:

dict[str, Any]

Returns:

Dictionary of parameter names and values

initialize_from_flow_curve(gamma_dot, sigma)

Initialize parameters from flow curve data.

Smart initialization strategy: 1. τ_y0 from low-rate plateau (stress intercept) 2. K_HB, n_HB from high-rate power-law region 3. G estimated from stress overshoot if available 4. Fluidity parameters from viscosity scaling

Parameters:

  • gamma_dot (ndarray) – Shear rate data (1/s)

  • sigma (ndarray) – Shear stress data (Pa)

Return type:

None

property pcov_

Parameter covariance matrix from NLSQ fit.

Returns:

ndarray of shape (n_params, n_params), or None if not fitted

property popt_

Optimal parameter values from NLSQ fit.

Returns:

ndarray of shape (n_params,), or None if not fitted

precompile(test_mode='relaxation', X=None, y=None)

Precompile NLSQ residual functions to eliminate JIT cold-start.

Triggers JIT compilation by running a minimal fit (max_iter=1) with dummy data. The model parameters are reset afterwards so the model is left in its original state.

This is useful for interactive sessions or benchmarks where the ~870ms first-fit JIT overhead should be excluded.

Parameters:

  • test_mode (str) – Test mode to precompile for (default: ‘relaxation’).

  • X (Optional[numpy.typing.ArrayLike]) – Optional input data for shape inference. If None, uses a 10-point logspace array.

  • y (Optional[numpy.typing.ArrayLike]) – Optional output data. If None, generates ones matching X.

Return type:

float

Returns:

Compilation time in seconds.

Example

>>> model = Maxwell()
>>> t = model.precompile(test_mode='relaxation')
>>> print(f"Compiled in {t:.2f}s")
>>> model.fit(X, y)  # No JIT overhead
precompile_bayesian(X=None, y=None, test_mode=None, num_chains=4)

Precompile NUTS kernel to eliminate JIT overhead in subsequent calls.

Triggers JIT compilation of the NumPyro model by running a minimal sampling (1 warmup, 1 sample). This caches the compiled kernel so that subsequent fit_bayesian() calls are 2-5x faster.

Parameters:

  • X (ndarray | RheoData | None) – Sample input data for determining array shapes. If None, uses a default 10-point linspace [0.01, 100].

  • y (ndarray | None) – Sample output data. If None, generates dummy data.

  • test_mode (str | TestModeEnum | None) – Test mode to precompile for. If None, defaults to ‘relaxation’.

Returns:

Compilation time in seconds.

Return type:

float

Example

>>> model = Maxwell()
>>> compile_time = model.precompile_bayesian(test_mode='relaxation')
>>> print(f"Compiled in {compile_time:.1f}s")
>>> # Now fit_bayesian() will be faster
>>> result = model.fit_bayesian(X, y)  # No compilation overhead
predict(X, test_mode=None, deformation_mode=None, poisson_ratio=None, **kwargs)

Make predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • test_mode (str | None) – Optional test mode (‘oscillation’, ‘relaxation’, ‘creep’, ‘flow’). If provided, sets model’s test_mode before prediction. Useful for data generation without fitting.

  • deformation_mode (str | DeformationMode | None) – Optional deformation mode for output conversion. If None, uses the mode stored from fit(). If tensile, converts G* predictions to E* space.

  • poisson_ratio (float | None) – Poisson’s ratio for conversion. If None, uses value stored from fit() (default 0.5).

  • **kwargs – Additional arguments passed to the internal _predict method.

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions (in E* space if deformation_mode is tensile)

predict_normal_stresses(gamma_dot)

Predict first and second normal stress differences.

The Saramito model predicts non-zero N₁ from its upper-convected Maxwell foundation.

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)

Return type:

tuple[ndarray, ndarray]

Returns:

  • N1 (np.ndarray) – First normal stress difference τ_xx - τ_yy (Pa)

  • N2 (np.ndarray) – Second normal stress difference τ_yy - τ_zz (Pa) Note: N2 = 0 for simple UCM models

property relaxation_time: float

Get characteristic relaxation time λ = 1/(G*f_age).

This is the relaxation time in the aged (rest) state.

Returns:

Relaxation time (s)

Return type:

float

sample_prior(num_samples=1000, seed=None)

Sample from prior distributions over parameter bounds.

Samples from uniform prior distributions defined by parameter bounds. This is useful for prior predictive checks and understanding the prior’s influence on the posterior.

Parameters:

  • num_samples (int) – Number of samples to draw from prior (default: 1000)

  • seed (int | None) – Random seed for reproducibility (default: None)

Return type:

dict[str, ndarray]

Returns:

Dictionary mapping parameter names to arrays of prior samples. Each array has shape (num_samples,) and dtype float64.

Raises:

Example

>>> model = MyModel()
>>> prior_samples = model.sample_prior(num_samples=500, seed=42)
>>> print(prior_samples["a"].shape)  # (500,)
score(X, y)

Compute model score (R² by default).

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – True target values

Return type:

float

Returns:

Model score (R² coefficient)

set_params(**params)

Set model parameters.

Parameters:

**params – Parameter names and values

Return type:

BaseModel

Returns:

self for method chaining

to_dict()

Serialize model to dictionary.

Return type:

dict[str, Any]

Returns:

Dictionary representation of model

property weissenberg_number: float | None

Get Weissenberg number Wi = λ * γ̇ (for steady/startup tests).

Returns:

Weissenberg number if gamma_dot is set, else None

Return type:

float or None

parameters: ParameterSet