In-memory Python solver

The simcoon.solver package drives the C++ material-point solver directly from Python: the loading path is defined with Block and StepMeca / StepThermomeca objects, and the results come back as numpy arrays — no path.txt, output.dat or result files involved. Since simcoon 2.0 this package is sim.solver, JSON is the only file format it reads, and the pre-2.0 text inputs are converted once with scripts/legacy_to_json.py.

Quick start

A uniaxial tension test on an elastic isotropic material:

import numpy as np
from simcoon import solver

step = solver.StepMeca(
    control=['strain'] + ['stress'] * 5,   # E11 driven, lateral stress-free
    value=[0.01, 0, 0, 0, 0, 0],           # targets, Voigt order [11,22,33,12,13,23]
    time=1.0, ninc=100,
)
res = solver.solve(step, "ELISO", [70000., 0.3, 1.E-5], nstatev=1)

stress = res["Stress"]     # Cauchy stress history, shape (6, N)
strain = res["Strain"]     # strain history (logarithmic for the default rate), shape (6, N)

Results follow the fedoo DataSet conventions — components first, one column per increment — so they interoperate directly with fedoo utilities (e.g. fedoo.util.voigt_tensors.StressTensorList(res["Stress"])). Available fields include Stress (Cauchy), Kirchhoff, PKII, Strain (the strain of the objective rate: logarithmic for the default rate, Almansi for "truesdell"), LogStrain (\(\ln\mathbf{V}\) from \(\mathbf{F}\), exact whatever the rate), GreenLagrange, F, R, DR ((3, 3, N)), TangentMatrix ((6, 6, N)), Statev, Wm, Time, Temp and, for thermomechanical runs, Q, r, Wt and the coupled tangents dSdE, dSdT, drdE, drdT.

Loading control

control sets each component to 'strain' (kinematically driven) or 'stress' (statically driven); mixed control is solved by Newton-Raphson.
Block(control_type=...) selects the strain/stress measures: 'small_strain' (default), 'green_lagrange' (PKII control), 'logarithmic' (Kirchhoff control), 'biot', or the fully kinematic 'F' / 'gradU' (9 components of the deformation gradient).
solve(corate=...) selects the objective rate for the finite-strain control types: 'jaumann', 'green_naghdi', 'logarithmic', 'logarithmic_R' (default), 'truesdell', 'logarithmic_F'. 'logarithmic_R' transports by the exact polar rotation increment \(\Delta\mathbf{R} = \mathbf{R}_1\mathbf{R}_0^T\), for which the tangent transport is exact — including with rotated internal-variable history (plasticity at finite rotation).
mode='sinusoidal' interpolates the step sinusoidally instead of linearly; mode='tabular' follows a user table passed in memory:

  t = np.linspace(0.01, 1.0, 100)
  e11 = 0.015 * np.sin(np.pi * t)
  step = solver.StepMeca(control=['strain'] + ['zero'] * 5,
                         mode='tabular',
                         tabular=np.column_stack([t, e11]))

In a saved path the table is the one input that is not JSON: ``save_path_json``
writes it as ``<stem>_tab<k>.csv`` next to the JSON (``#`` header naming the
columns, one row per increment) and the step's ``"tabular"`` entry holds that
filename; ``load_path_json`` reads it back, comma- or whitespace-separated.

Cyclic loading repeats the steps of a block: Block(steps=[...], ncycle=10). Tabular steps cannot be cycled (their time column is absolute); unroll the cycles into explicit steps instead.
A rotation rate can be superimposed on the mixed finite-strain control types through StepMeca(BC_w=...) (3x3 spin matrix).

Thermomechanical loading

StepThermomeca activates the coupled heat equation (block type 2, small strain), with thermal_control set to 'temperature' (ramp to T_final), 'heat_flux' (prescribed Q) or 'convection' (0D convection with coefficient q_conv):

step = solver.StepThermomeca(control=['stress'] * 6, value=[0.] * 6,
                             T_final=340., ninc=50)
# ELISO thermomechanical props: rho, c_p, E, nu, alpha
res = solver.solve(step, "ELISO", [1.E-9, 1., 70000., 0.3, 1.E-5], nstatev=1)

Mechanical steps (StepMeca) also accept T_final: the temperature then ramps as an imposed condition of the mechanical problem (thermal expansion without the heat equation).

Tabular thermomechanical steps support the three thermal controls: with thermal_control='temperature' the table carries a T column when tabular_T=True (constant temperature otherwise); with 'heat_flux' the thermal column is the prescribed flux Q; with 'convection' there is no thermal column and q_conv applies.

For finite-element couplers, the point-wise thermomechanical UMAT batch entry sim.umat_T(...) complements sim.umat(...); it returns (sigma, statev, Wm, Wt, r, dSdE, dSdT, drdE, drdT).

Mean-field composites

The mean-field models (MIHEN, MIMTN, MISCN, MIPLN) take their sub-phases in memory, as Ellipsoid or Layer objects passed to solve(phases=...) or sim.L_eff(..., phases=...) (see Use the solver, phases); their own props carry only the scheme’s settings, [mp, np, n_matrix] for Mori-Tanaka (see Constitutive model (UMAT) catalog). Every orientation in those objects is a simcoon.Rotation: the material frame of a phase (material_orientation) and the geometry of an inclusion or a layer (geometry_orientation). Any of these forms is accepted and coerced:

from simcoon.solver.micromechanics import Ellipsoid, as_rotation, euler_angles

fibre = Ellipsoid(umat_name="ELISO", concentration=0.2, nstatev=1,
                  props=[50000., 0.3, 0.], a1=50.,
                  geometry_orientation=sim.Rotation.from_rotvec([0, 0, np.pi / 4]))
fibre.geometry_orientation = (45., 0., 0.)                       # Euler angles, degrees
fibre.geometry_orientation = {"psi": 45., "theta": 0., "phi": 0.}  # the JSON form
euler_angles(fibre.geometry_orientation)   # {'psi': 45.0, 'theta': 0.0, 'phi': 0.0}

The Euler angles are the 'zxz' sequence the C++ side reads, in degrees: as_rotation((psi, theta, phi)) is Rotation.from_euler('zxz', [psi, theta, phi], degrees=True) (scipy’s extrinsic zxz), applied actively. A phase at that orientation responds with R.apply_stiffness(L_local), which the test suite pins against the solver and L_eff. The JSON files and the dicts handed to the extension keep the angles; euler_angles writes them back, with the usual caveat that the decomposition is not unique at theta = 0 (the z angles merge into psi). solve(orientation=...) and sim.L_eff(umat_name, props, nstatev, orientation=..., phases=...) take the same forms for the frame of the material or of the whole RVE; L_eff also takes the phase objects directly.

An orientation distribution splits one phase into phases rotated about a direction, with concentrations following the ODF:

from simcoon.solver.micromechanics import Peak, discretize_odf

peaks = [Peak(method=3, mean=90., s_dev=10.)]     # Gaussian, degrees
phases = discretize_odf([matrix, fibre], num_phase=1, peaks=peaks, nphases=18,
                        axis=(0., 0., 1.), angle_range=(0., 180.))
L = sim.L_eff("MIMTN", props, nstatev, phases=phases)

The k-th copy is rotated by alpha_k about axis on top of the orientations the phase already has (Rotation.from_rotvec(alpha_k * axis) * orientation), the material frame following unless rotate_material=False; alpha_k runs from angle_min in nphases equal steps and each copy takes the Simpson integral of the density over its step, normalised to the parent’s concentration. The density is a director distribution, periodic over 180 degrees: the default half turn is right when a half turn about axis maps the inclusion onto itself (axis along or normal to a principal axis of the inclusion); about any other direction give angle_range=(0., 360.). The pre-2.0 Euler sweeps are the cases axis=(0, 0, 1) (psi or phi) and axis=(1, 0, 0) (theta) from an unrotated phase. The density itself is sim.get_densities_ODF(x, peaks), the sum of the Peak profiles (method: 1 standard-deviation kernel, 2 hard cut-off, 3 Gaussian, 4 Lorentzian, 5 pseudo-Voigt, 6 Pearson VII, 7 uniform; formulas in the class), useful to plot a distribution before discretising it; Peak.density(x) without periodic is the plain profile, for a distribution of a parameter rather than of an orientation.

JSON configuration

Materials and loading paths round-trip through JSON (save_material_json(), save_path_json(), load_simulation_json()):

solver.save_material_json("material.json", "ELISO", [70000., 0.3, 1.E-5], 1)
solver.save_path_json("path.json", [solver.Block(steps=[step])], T_init=293.15)
res = solver.solve(**solver.load_simulation_json("material.json", "path.json"))

Results can be persisted with res.save("run.npz") / SolverResults.load("run.npz"), or flattened with res.to_dataframe().

Legacy text inputs

JSON is the only format simcoon reads or writes since 2.0: nothing in the package parses path.txt, material.dat, tab_file_<n>.txt or N<kind><n>.dat. A pre-2.0 data directory is converted once with the migration script shipped in the repository (not in the package), the tables its mode-3 steps referenced being rewritten as <stem>_tab<k>.csv files next to the path JSON:

python scripts/legacy_to_json.py data   # writes path.json, material.json, ellipsoids<N>.json, ...

res = solver.solve(**solver.load_simulation_json("data/material.json", "data/path.json"))

Solver parameters

solve() exposes the numeric controls of the adaptive Newton loop as keyword arguments: precision (default 1e-6), maxiter/miniter, div_tnew_dt/mul_tnew_dt (time-step cut/growth factors), inforce, lambda_solver (penalty stiffness of strain-driven components), plus tangent_mode ('none', 'continuum', 'algorithmic' — default) and solver_type.

Constitutive laws written in Python

umat_name also accepts a law object implementing simcoon.PythonUMAT (numpy, PyTorch, …): it is registered under the PYEXT name for the duration of the call and integrated by the C++ solver like a built-in kernel; props and nstatev default to the object’s attributes. See Constitutive laws in Python (PYEXT).

API reference

class simcoon.solver.StepMeca(control: str | Sequence[str | int] = 'strain', value: Sequence[float] | None = None, time: float = 1.0, ninc: int = 100, mode: str | int = 'linear', Dn_init: float = 1.0, Dn_mini: float = 0.001, BC_w: Sequence[Sequence[float]] | None = None, T_final: float | None = None, tabular: ndarray | None = None, tabular_T: bool = False)

A mechanical loading step.

Parameters:

control (str or sequence of str) – Per-component control: ‘strain’ (kinematic) or ‘stress’ (static); ‘zero’ holds a component at zero for tabular steps. A single string applies to all components.
value (array-like, optional) – Absolute target values at the end of the step, per component, in Voigt order (6 components; 9 row-major for control types ‘F’/’gradU’). Not used for tabular steps.
time (float) – Duration of the step.
ninc (int) – Number of increments (linear and sinusoidal modes).
mode (str or int) – ‘linear’, ‘sinusoidal’ or ‘tabular’.
Dn_init (float, optional) – Initial and minimal sub-increment fraction of one increment (adaptive stepping). Defaults: 1.0 and 1e-3.
Dn_mini (float, optional) – Initial and minimal sub-increment fraction of one increment (adaptive stepping). Defaults: 1.0 and 1e-3.
BC_w (array-like, optional) – 3x3 spin matrix (rotation rate) applied during the step, for the mixed finite-strain control types (‘green_lagrange’, ‘logarithmic’, ‘biot’).
T_final (float, optional) – Target temperature at the end of the step (temperature ramp applied to the mechanical UMAT). None (default) holds the temperature.
tabular (ndarray, optional) – Mode-3 table, one row per increment. Columns: [time, (thermal column, see tabular_T and StepThermomeca.thermal_control), controlled components in Voigt order]. Required when mode=’tabular’. The time column is ABSOLUTE simulation time and must continue from the previous step’s end time (a table restarting at 0 after an earlier step is rejected: it would produce negative time increments). Tabular steps cannot be cycled (Block.ncycle must be 1).
tabular_T (bool) – Whether the tabular table contains a temperature column (after the time column). Default False (constant temperature). For thermomechanical heat-flux steps the thermal column is the flux Q and is always required, regardless of this flag.

T_end(T_hold: float) → float: Temperature at the end of the step (for chaining T_final=None steps).

to_dict(control_type: int, T_hold: float) → dict

Marshal to the dict consumed by simcoon._core.solver_run.

Parameters:

control_type (int) – The control type of the enclosing block (drives the 6/9 sizing).
T_hold (float) – Temperature to hold when T_final is None (running block value).

class simcoon.solver.StepThermomeca(control: str | Sequence[str | int] = 'strain', value: Sequence[float] | None = None, time: float = 1.0, ninc: int = 100, mode: str | int = 'linear', Dn_init: float = 1.0, Dn_mini: float = 0.001, BC_w: Sequence[Sequence[float]] | None = None, T_final: float | None = None, tabular: ndarray | None = None, tabular_T: bool = False, thermal_control: str = 'temperature', Q: float = 0.0, q_conv: float = 0.0)

A thermomechanical loading step (coupled heat equation).

In addition to the mechanical control of StepMeca:

Parameters:

thermal_control (str) – ‘temperature’ (ramp to T_final), ‘heat_flux’ (prescribed flux Q) or ‘convection’ (0D convection Q = -q_conv (T - T_init)).
Q (float) – Prescribed heat flux (thermal_control=’heat_flux’).
q_conv (float) – Convection coefficient rho*c_p/tau (thermal_control=’convection’).

T_end(T_hold: float) → float: Temperature at the end of the step; flux/convection steps leave the chained hold temperature unchanged (the reached T is solution-dependent).

class simcoon.solver.Block(steps: ~typing.List[~simcoon.solver.blocks.StepMeca] = <factory>, control_type: str | int = 'small_strain', ncycle: int = 1)

A loading block: a sequence of steps repeated ncycle times.

Parameters:

steps (list of StepMeca / StepThermomeca) – The loading steps. If any step is a StepThermomeca, the block is thermomechanical (coupled heat equation) and all steps must be.
control_type (str or int) – Loading control (see CONTROL_TYPES). Thermomechanical blocks only support ‘small_strain’. A run is either all small strain or all finite strain: a finite-strain block restarts from F = I after a small-strain one.
ncycle (int) – Number of repetitions of the step sequence.

to_dict(T_hold: float) → dict: Marshal to the dict consumed by simcoon._core.solver_run.

Solve a homogeneous loading path with the C++ simcoon solver, in memory.

Parameters:

blocks (Block, StepMeca or sequence of them) – The loading path. Bare steps are wrapped in a small-strain Block.
umat_name (str, PythonUMAT or callable) – Constitutive model: either the name of a built-in model (5 characters, e.g. ‘ELISO’, ‘EPICP’, ‘MODUL’) or a constitutive law written in Python (a simcoon.PythonUMAT instance, or any callable with its integrate keyword signature). A Python law is registered under the PYEXT name for the duration of the call and integrated by the C++ solver exactly like a built-in kernel.
props (array-like, optional) – Material properties. Required for a built-in model; defaults to the props attribute of a Python law.
nstatev (int, optional) – Number of internal state variables. Required for a built-in model; defaults to the nstatev attribute of a Python law.
T_init (float) – Initial temperature.
corate (str, int or None) – Objective rate for the finite-strain control types (see CORATE_TYPES). Default (None): ‘logarithmic_R’ — the exact polar rotation, whose frame increment DR = R1 R0^T makes the tangent transport exact even with rotated internal-variable history (the XBM ‘logarithmic’ rate keeps a small tangent residual there). MODUL additionally requires log_R under NLGEOM (the modular Hencky composition).
tangent_mode (str or int) – Tangent operator mode: ‘none’, ‘continuum’, ‘algorithmic’ (default) or ‘closest_point’ (reserved).
solver_type (int) – 0 = classic Newton-Raphson (default), 1 = RNL (control_type 1 only).
orientation (simcoon.Rotation, dict or sequence of 3 floats) – Orientation of the material frame: a Rotation, or its Euler angles (psi, theta, phi) in degrees (see as_rotation()); applied actively, material frame to global frame.
phases (sequence, optional) – Sub-phases of a mean-field model (MIMTN, MISCN, MIHEN, MIPLN): the Ellipsoid / Layer objects of simcoon.solver.micromechanics, or the dicts they convert to. Leave it None for every single-phase model.
record_tangent (bool) – Capture the tangent operator history (‘TangentMatrix’ or the coupled thermomechanical tangents).
raise_on_abort (bool) – Raise a RuntimeError when the solver aborts early (status != 0) instead of returning the partial history. The solver aborts when the Newton loop does not converge at the minimal increment, or when the increment falls below Dn_mini with inforce=0.
**params – Numeric solver controls forwarded to the C++ loop: div_tnew_dt, mul_tnew_dt, miniter, maxiter, inforce, precision, lambda_solver (penalty stiffness of the strain-driven components).

Returns:

History of the converged increments (fedoo-style data layout).

Return type:

SolverResults

class simcoon.solver.SolverResults(raw: Dict[str, ndarray], finite_blocks=None)

Dict-like access to the solver history.

scalar_data

Scalar histories, shape (N,): ‘Time’, ‘Temp’, ‘Block’, ‘Cycle’, ‘Step’, ‘Inc’; thermomechanical runs add ‘Q’ (heat flux) and ‘r’ (heat source).

Type:: dict

field_data

Tensor histories, components-first: ‘Stress’ (Cauchy, (6, N)), ‘Kirchhoff’, ‘PKII’, ‘Strain’ (the strain integrated with the objective rate, (6, N): ln V for the logarithmic rates, the Almansi strain for ‘truesdell’), ‘LogStrain’ (ln V computed from F whatever the rate; the small strain on small-strain blocks), ‘GreenLagrange’ ((6, N), computed from F), ‘Statev’ ((nstatev, N)), ‘Wm’ ((4, N)), ‘F’, ‘R’, ‘DR’ ((3, 3, N)); ‘TangentMatrix’ ((6, 6, N)) for mechanical runs; thermomechanical runs add ‘Wt’ ((3, N)) and the coupled tangents ‘dSdE’ ((6, 6, N)), ‘dSdT’ ((6, N)), ‘drdE’ ((6, N)), ‘drdT’ ((N,)).

Type:: dict

status

0 if the simulation ran to completion, 1 on early abort (the recorded history is then partial).

Type:: int

get_data(key: str) → ndarray: fedoo-style accessor (alias of __getitem__).

classmethod load(filename: str) → SolverResults: Load a SolverResults previously written by save().

save(filename: str) → None: Save all histories to a compressed npz archive.

to_dataframe(): Flatten scalar and 6-component histories to a pandas DataFrame.