Physics informed neural networks for continuum micromechanics

In plain words

Engineers want to know how a solid deforms under load and where the stress piles up. A physics-informed neural network learns that answer by trying to satisfy the governing equations at many sample points, without measured data. This paper applies the idea to materials made of a stiff phase inside a soft one, including a slice of a scanned wood-plastic composite. Sharp borders between the two phases are hard for a single network to fit, because one network has to describe the whole piece at once.

The authors test two remedies. One puts more sample points where the error is large. The other cuts the piece into patches and gives each patch its own small network, with rules that make neighbouring patches agree on their borders. The patch approach works better. The error-driven sampling does not converge, and the authors say it is not yet a faster solver than established methods for forward problems.

TL;DR

The network outputs two displacements and three stress components, so no second derivatives are needed, and the boundary conditions are built into the output. Training uses full BFGS on the local equation residuals plus one global internal-versus-external work term. On a square plate with one stiff circular inclusion, adaptive collocation fails to converge, while a conservative PINN (cPINN) split into a 4x4 grid gives the highest accuracy with regular sampling §7. The released code is GPL-3.0 and differs from the paper in several settings. All results are author-reported, and the micro-CT experiment cannot be rerun from the release.

Problem

The governing equations are the static, linear elastic balance of momentum without body forces, with symmetric stress, small-strain kinematics and a linear constitutive law §2 eq. 5 eq. 7 eq. 8. The domain is a zero-centred square eq. 1. The paper assumes plane strain fig. 3.

A PINN uses one global ansatz for the whole domain, while FEM uses local ansatz functions §3. For material jumps this has three costs. Sharp phase transitions can make gradients explode during optimisation. A displacement-only formulation needs second derivatives of the displacement §4. Free stress boundaries are hard to impose with the hard-constraint trick §4. A single network also shows its largest errors at the phase transitions §5.1. For real micro-CT microstructures no analytic material law exists, so the material field must be learned §6.1.

Method

Network. A densely connected feed-forward network with a constant number of units per layer and a linear last layer §3.1 eq. 14 eq. 15. The paper uses the Swish activation with a slope parameter β eq. 16. The value of β is not stated. Biases are absorbed into the weight vectors §3.1.

Outputs. The network returns five values: u_x, u_y, σ_xx, σ_yy and σ_xy §4 eq. 25. Stress is an output, so the balance law is evaluated with first derivatives only. The symmetry of the stress tensor removes σ_yx §2 eq. 6.

Loss. The PINN is trained by collocation: the loss is the mean squared residual over the domain points, plus a boundary term in the soft-constraint case §3.2.1 eq. 19 eq. 20 eq. 21. Spatial derivatives come from automatic differentiation §3.2.1.

Hard boundary conditions. Each field is written as a prescribed part plus a smooth distance function times a raw network output §3.2.2 eq. 22. With hard constraints the boundary term disappears and the loss keeps only the residuals eq. 23 §3.2.2. The distance factors for the plate are given in §4.3.1; the σ_xy factor is discussed under Paper vs code.

Residual loss for linear elasticity. The local loss has five mean-squared terms (two balance equations, three constitutive equations) and one global work term §4.1. The work term is the squared difference between the internal and external work over the whole domain, with the factors shown in §4.1. The paper uses no weights on the local terms §4.1, remark.

Scaling. Each quantity is divided by a reference constant before the loss is computed, and predictions are multiplied back afterwards §4.2. The paper says the constants are the maxima of each quantity over the domain. The scaling changes only the loss balance, not the PDE §4.2.

Optimiser. Full BFGS, not L-BFGS and not Adam §4.1, remark §4.3.1.

Adaptive collocation (algorithm 1). Starting from a fine regular grid, each cycle keeps a coarse regular grid, draws uniform random candidates, computes the local error at each candidate without the work term, and keeps the candidates with the largest error §5.1 alg. 1. The regular and adaptive sets are then combined, and the ratio γ = n_reg / n_ada sets how many points come from each §5.1. Gradients are clipped to a threshold α during training §5.1. The paper does not state α.

Domain decomposition (cPINN). The domain is split into subdomains, each with its own network §5.2. Neighbouring subdomains must agree on the displacement and on the traction across their shared edge §5.2. The paper adds these agreement conditions as squared residuals with a constant weight ψ §5.2. The hard boundary constraints must be adapted for subnetworks that share part of the outer boundary §5.2.

Material field. For the single inclusion, λ and μ follow a tanh profile of the radius, with constants c1, c2, c3 and a steepness δ that the paper does not give §4.3.2. For a real micro-CT slice, a separate network is trained on the binarised microstructure, with a tanh-scaled output that keeps λ and μ inside their physical range §6.1. Its training is described as lasting minutes §6.1.

# Hard-constrained PINN outputs for the plate  [§4] [§4.3.1]
# Input: coords (N, 2) in [-1, 1]^2; output: five fields of shape (N,)
def pinn_fields(coords, sigma_bar):
    x = coords[:, 0]
    y = coords[:, 1]
    raw = network(coords)                       # (N, 5)  [§3.1]
    # Prescribed part plus distance factor times raw output  [§3.2.2]
    u_x = (-1 - x) * raw[:, 0]                  # zero on x = -1
    u_y = (-1 - y) * raw[:, 1]                  # zero on y = -1
    sigma_xx = sigma_bar + (1 - x) * raw[:, 2]  # sigma_bar on x = 1
    sigma_yy = (1 - y) * raw[:, 3]              # zero on y = 1
    # Paper's factor; see Paper vs code for the code's variant  [§4.3.1]
    sigma_xy = (1 - x)**2 * (1 - y)**2 * raw[:, 4]
    return u_x, u_y, sigma_xx, sigma_yy, sigma_xy

# Local PINN loss with the global work balance  [§4.1]
# lam, mu: Lame fields at the points  (N,)
def pinn_loss(coords, lam, mu, sigma_bar):
    u_x, u_y, sigma_xx, sigma_yy, sigma_xy = pinn_fields(coords, sigma_bar)
    # Balance of linear momentum, no body force  [eq. 5]
    r_x = grad_x(sigma_xx) + grad_y(sigma_xy)
    r_y = grad_x(sigma_xy) + grad_y(sigma_yy)
    # Constitutive law with the small-strain kinematics  [eq. 7] [eq. 8]
    trace = grad_x(u_x) + grad_y(u_y)
    c_xx = lam * trace + 2 * mu * grad_x(u_x) - sigma_xx
    c_yy = lam * trace + 2 * mu * grad_y(u_y) - sigma_yy
    c_xy = mu * (grad_y(u_x) + grad_x(u_y)) - sigma_xy
    # Mean squared residuals over the n_d domain points  [§4.1]
    loss = mean(r_x**2) + mean(r_y**2)
    loss = loss + mean(c_xx**2) + mean(c_yy**2) + mean(c_xy**2)
    # Internal work over the domain; L is the edge length  [§4.1]
    # integrand: sxx*dux/dx + syy*duy/dy + sxy*(dux/dy + duy/dx)
    work_int = (L**2 / (2 * n_d)) * sum(internal_work_integrand)
    # External work on the traction edge, 1/n_b per boundary point  [§4.1]
    work_ext = (1 / n_b) * sum(sigma_xx * u_x + sigma_xy * u_y)
    loss = loss + (work_int - work_ext)**2
    return loss

# Adaptive collocation training  [alg. 1] [§5.1]
def adaptive_training(net, n_fine, n_iter, n_reg, n_ada, alpha):
    # Pre-training on the fine regular grid  [alg. 1]
    for step in range(n_fine):
        train(net, regular_grid())
    for cycle in range(n_iter):
        # Coarse regular grid, held fixed across cycles  [§5.1]
        grid_points = coarse_regular_grid(n_reg)
        # Uniform random candidates in the square  [§5.1]
        candidates = uniform(-L / 2, L / 2, size=(n_rand, 2))
        # Local error at each candidate, without the work term  [§5.1]
        errors = pde_error(net, candidates)
        # Keep the n_ada candidates with the largest error  [§5.1]
        adaptive_points = candidates[argsort(-errors)[:n_ada]]
        points = concat(grid_points, adaptive_points)
        # BFGS with gradients clipped to norm alpha  [§5.1]
        train(net, points, clip_norm=alpha)
    return net

# Conservative PINN loss: one subnetwork per subdomain  [§5.2]
# Shared edges counted once; the paper leaves this open  [§5.2]
# psi = 20 in the experiments  [§5.3]
def cpinn_loss(subnets, edges, psi):
    total = 0
    for net in subnets:
        total = total + pinn_loss(net)          # local loss  [§4.1]
    for edge in edges:
        # Horizontal edge between I and II: u_x and both stress
        # components sigma_xx, sigma_xy must agree  [§5.2]
        if edge.is_horizontal:
            total = total + psi * mean((u_x[I] - u_x[II])**2)
            total = total + psi * mean((sigma_xx[I] - sigma_xx[II])**2)
            total = total + psi * mean((sigma_xy[I] - sigma_xy[II])**2)
        # Vertical edge: u_y, sigma_yy, sigma_xy  [§5.2]
        else:
            total = total + psi * mean((u_y[I] - u_y[II])**2)
            total = total + psi * mean((sigma_yy[I] - sigma_yy[II])**2)
            total = total + psi * mean((sigma_xy[I] - sigma_xy[II])**2)
    return total

Hyperparameters

Results

Author-reported. Setup: 2D square [-1, 1]^2 in linear elasticity, with the plate and the single-inclusion experiments of §4.3, the convergence study of §5.3, and a micro-CT slice of a wood-plastic composite of §6. The residual R is the sum of absolute scaled residuals §4.3.1. These values are scaled, so they are not directly physical errors.

Mean R against the number of collocation points n_d, single inclusion fig. 8. Values are the plotted data:

The text says the cPINN has lower error than the PINN throughout §5.3. The table agrees. The adaptive variants do not converge from 64^2 upwards §5.3. The text also says the adaptive PINN is lower than the PINN in the first half of the study. The table shows this only at 16^2 and 32^2; at 8^2 the adaptive PINN is higher fig. 8.

Paper vs code

The official code is at ahenkes1/HENKES_PINN, pinned at the commit in the frontmatter. The owner handle matches the first author's name, and the repository README and CITATION file name the authors README.md. The repository is GPL-3.0 (copyleft). This card cites paths and lines only and quotes no code. The arXiv comment points to the journal version. The paper text used here is arXiv v2.

• Activation. Paper: Swish, β not stated §3.1 eq. 16. Code: tanh main.py.
• Hidden layers. Paper: n_L = 4 §4.3.1, and Figure 2 draws three hidden layers fig. 2. Code: four hidden Dense layers of the chosen width, then a 5-unit output layer cpinn.py main.py. By our count from these layer sizes, four hidden layers of 64 units give about 13.0 x 10^3 weights, matching the paper's figure §5.3. Three hidden layers would give about 8.8 x 10^3. The code matches the count, not the drawing.
• cPINN width. Paper: width chosen to keep about 13 x 10^3 weights in total §5.3. Code: width = 64 / split, so 16 units for 4x4 main.py. By our count this gives about 15.2 x 10^3 weights in total.
• CT cPINN width. Paper: n_u = 90 §6.2.2. Code: width is fixed by the formula above; there is no flag for it main.py.
• Collocation. Paper: n_d = 128^2 §4.3.1. Code: default --points 32 per dimension main.py; final prediction uses 128 per dimension main.py. The regular grid includes both edges of each subdomain dataset.py.
• Stopping. Paper: plate stops after 138 iterations on gradient size §4.3.1. Code: the non-adaptive BFGS call sets all tolerances to zero, and the iteration count is --iter, default 1000 main.py [code: src/main.py#L710]. Line search is capped at 50 steps [code: src/main.py#L307].
• Gradient clipping. Paper: used in the adaptive scheme only, threshold not stated §5.1. Code: the gradient passed to BFGS is clipped to norm 1 in every run, adaptive or not bfgs.py. The loss values used by the line search are not clipped bfgs.py.
• Shear-stress factor. Paper: (1 − x)^2 (1 − y)^2 §4.3.1. This leaves σ_xy nonzero on x = −1 and y = −1, although the boundary conditions in Figure 3 set σ_xy = 0 on all four edges fig. 3. Code: (x^2 − 1)(y^2 − 1), which vanishes on all four edges cpinn.py. The code follows the figure.
• σ_yy factor. Paper: (1 − y) §4.3.1. Code: (−1 + y), the same zero set with the sign absorbed into the output cpinn.py.
• Scaling. Paper: each quantity divided by its maximum over the domain §4.2. Code: fixed constants σ_c = 0.025, u_c = L = 1, λ_c = 14285.7, μ_c = 3571.4 cpinn.py. Predictions are multiplied back by σ_c / λ_c for displacements and by σ_c for stresses cpinn.py [code: src/cpinn.py#L744-L748].
• Lamé constants. Paper: derived from E and ν §4.3.1 §4.3.2, without the conversion printed. Code: the stiff and soft values are hard-coded and equal λ = Eν / ((1 + ν)(1 − 2ν)) and μ = E / (2(1 + ν)) for E = 10^4 and 1.5 x 10^3 with ν = 0.4 cpinn.py. The same relations appear in the dataset helper dataset.py.
• Inclusion profile. Paper: c1, c2, c3, δ not stated §4.3.2. Code: radius 0.4 and δ = 0.01 cpinn.py. The code's profile is always used. There is no switch to the homogeneous plate of §4.3.1.
• Work term. Paper: squared difference of integrals with L^2/(2 n_d) and 1/n_b factors §4.1. Code: relative squared error, divided by the larger of the two magnitudes, with weight 1 cpinn.py [code: src/cpinn.py#L528].
• Loss weights. Paper: none for the local PINN terms §4.1, remark; ψ = 20 for cPINN interface terms §5.3. Code: 1 for local terms [code: src/cpinn.py#L526-L528] and 20 for all interface and flux terms [code: src/cpinn.py#L530] cpinn.py. These agree.
• Interface pairs. Paper: each subnetwork's loss includes the interface terms, and the total is a sum over subnetworks §5.2. The paper does not say whether a shared edge is counted once or twice. Code: each shared edge is added once cpinn.py.
• Adaptive schedule. Paper: a fine-grid pass and one adaptive cycle, each with 5 x 10^3 BFGS iterations §5.3 alg. 1. Code: two passes, a regular grid pass and one adaptive pass, both with --iter iterations main.py main.py main.py. Ratio γ = 2.25 main.py. Candidate pool of 128^2 points per subnetwork main.py main.py. Error for selection excludes the work term cpinn.py.
• Material network. Paper: trained network for the CT slice, with a tanh output map for λ and μ §6.1. Code: a branch exists, but main.py never passes a material network main.py and the CT data is not in the release cpinn.py. The μ branch adds 1 to the network output, but the λ branch does not cpinn.py [code: src/cpinn.py#L361-L365]. The paper's map shifts both by the same amount §6.1.
• Optimiser object. Compiled with Adam, but training calls only BFGS. Adam has no effect on the results main.py.

Pitfalls

• Scaled residuals. R and the norms in this card are in scaled units §4.2 §4.3.1. They cannot be compared across experiments that use other constants.
• Scaling choice changes training. The paper says scaling affects only the loss balance §4.2. The code's fixed constants differ from the paper's stated rule cpinn.py. Expect different optimisation paths and different final residuals.
• Algorithm 1 loop bound. The inner training loop of alg. 1 is written with the bound n_ada, which looks like a typo. The text says each cycle runs 5 x 10^3 BFGS iterations §5.3. Implement the text's count.
• Clipping formula. The displayed clipping rule in §5.1 has a misplaced brace. Implement the usual clip to norm α.
• Clipped gradients in BFGS. Clipping the gradient while the line search uses the unclipped loss can make the quasi-Newton model inconsistent with the objective. The code does this in every run bfgs.py. Test this before relying on the paper's iteration counts.
• Early stop on GPUs. The README says BFGS can stop before the requested iteration count on GPUs without double-precision support, and that mixed precision loss scaling is a workaround README.md. The run is single precision §4.3 main.py.
• Edge points in training. The regular grid includes the subdomain edges, with no offset dataset.py. The residual is therefore enforced on the boundary as well. Check whether your implementation expects interior points only.
• Adaptive collocation is fragile. It does not converge from 64^2 points upwards and is sensitive to γ and to the cycle count §5.3 fig. 8. The paper concludes that the scheme does not improve accuracy in the form studied, and that many variants of it exist §5.3.
• More subnetworks is not always better. The 5x5 split has a higher mean residual than 4x4 fig. 10. Each extra subnetwork adds interface terms that make the loss harder to optimise §5.3.
• Percentage claim. The text puts the reduction of the 4x4 cPINN's mean residual, relative to the plain PINN, at one quarter §5.3. The plotted values give 1.22 x 10^-1 against 1.49 x 10^-1, about 18% lower fig. 10. Use the figure values.
• Shared-boundary hard constraints. Subnetworks that share part of the outer boundary need adapted hard constraints §5.2. The code builds one distance function per subnetwork and edge case cpinn.py. Check every case if you rewrite it.
• Work term is global. It is not a per-point loss, so it cannot be used to select adaptive points §5.1.
• Unstated values. Activation β, initialisation, gradient clipping threshold, inclusion constants c1, c2, c3, δ, and the candidate count in the adaptive scheme are not in the paper §3.1 §4.3.2 §5.1. The code fills some of them in, as listed above.
• Parameter count claim. The paper gives about 2.7 x 10^4 weights for the CT cPINN §6.2.2. By our count from the layer sizes, one 90-unit subnetwork has about 2.5 x 10^4 weights and 16 subnetworks have about 4 x 10^5. The paper does not say which count it means.
• Compute. The paper says forward problems are not competitive with FEM or FFT in time, and that large problems can take hours on a GPU §7. Inverse problems are the stated next step §7.

Lineage

• Builds on: the PINN framework of Raissi et al. (2019) §3.2; the hard boundary construction of Berg and Nyström (2018) §3.2.2; the conservative PINN of Jagtap, Kharazmi and Karniadakis (2020) §5.2; and an adaptive collocation approach by Wight and Zhao (2020), which the paper says is related but different §5.1.
• Compared against: standard PINN, cPINN, adaptive PINN and adaptive cPINN on the single inclusion §5.3.
• Related extensions named by the authors: XPINNs and parallel cPINNs §7.