{'fraction_non_alphanumeric': 0.033099082667238884, 'fraction_numerical': 0.010112118448581085, 'mean_word_length': 4.784794414274631, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Composable Energy Surrogates for PDE Order Reduction

Meta-materials are an important emerging class of engineered materials in which complex macroscopic behaviour--whether electromagnetic, thermal, or mechanical--arises from modular substructure. Simulation and optimization of these materials are computationally challenging, as rich substructures necessitate high-fidelity finite element meshes to solve the governing PDEs. To address this, we leverage parametric modular structure to learn component-level surrogates, enabling cheaper high-fidelity simulation. We use a neural network to model the stored potential energy in a component given boundary conditions. This yields a structured prediction task: macroscopic behavior is determined by the minimizer of the system's total potential energy, which can be approximated by composing these surrogate models. Composable energy surrogates thus permit simulation in the reduced basis of component boundaries. Costly ground-truth simulation of the full structure is avoided, as training data are generated by performing finite element analysis with individual components. Using dataset aggregation to choose training boundary conditions allows us to learn energy surrogates which produce accurate macroscopic behavior when composed, accelerating simulation of parametric meta-materials.


Introduction
Many physical, biological, and mathematical systems are most successfully modeled by partial differential equations (PDEs). Analytic solutions are rarely available for PDEs of practical importance; thus, computational methods to approximate PDE solutions are critical for tackling many problems in science and engineering. Finite element analysis (FEA) is one of the most widely used techniques for solving PDEs on spatial domains; the continuous problem is discretized and replaced by basis functions on a mesh.
The effectiveness of FEA and related methods is largely governed by the granularity of the discrete approximation, i.e., the fineness of the mesh. Complicated domains can require fine meshes that result in prohibitively expensive computations to solve the PDE. This problem is compounded when the task is one of parameter identification or design optimization. In these situations the PDE must be repeatedly solved in the inner loop of a bi-level optimization problem.
An important domain where this challenge is particularly relevant is in the modeling of mechanical meta-materials. Meta-materials are solids in which heterogeneous microstructure leads to rich spaces of macroscopic behavior. Meta-materials offer the possibility of achieving electromagnetic and/or mechanical properties that are otherwise impossible with homogenous materials and traditional design approaches (Poddubny et al., 2013;Cai & Shalaev, 2010;Bertoldi et al., 2017). In this paper we focus on the class of cellular mechanical meta-materials proposed by Overvelde & Bertoldi (2014), which promise new highperformance materials for soft robotics and other domains (see Sec 3). The simulation of these meta-materials is particularly challenging due to the complex cellular structure and the need to accurately capture small-scale nonlinear elastic behavior. Traditional finite element methods have limited ability to scale to these problems, and automated design of these promising materials demands accurate and efficient approximate solutions to the associated PDE.
We develop a framework which takes advantage of spatially local structure in large-scale optimization problems-here the minimization of energy as a function of meta-material displacements. Given that only a small subset of material displacements are of interest, we "collapse out" the remainder using a learned surrogate. Given a component with substructure defined by local parameters, the surrogate produces an accurate proxy energy in terms of the boundary alone. A single surrogate can be trained, then replicated and composed to predict the energy of a larger solid-which is the sum of energies over sub-components. This allows solving the PDE in the reduced basis of component boundaries by minimizing the sum of surrogate energies.
Other methods exist for amortizing the solution of PDEs.

arXiv:2005.06549v2 [cs.LG] 15 May 2020
Some of the most common approaches use neural networks to map from PDE parameters to solutions (Zhu et al., 2019;Nie et al., 2020) or construct reduced bases via solving eigenvalue problems or interpolating between snapshots (Berkooz et al., 1993;Chatterjee, 2000). These approaches typically require solving full systems to produce training data. Our framework uses the modular decomposition of energy to train surrogate models on data generated by querying the finite element "expert" on the energy in small components, avoiding performing FEA on large systems which are expensive to solve.
The remaining sections of this paper are as follows. Section 2 proposes a framework for collapsed-basis optimization when the objective decomposes as a sum of terms which each depend only on local variables. Section 3 provides a brief introduction to the cellular solids of interest and the PDEs governing their behavior. In Section 4, collapsedbasis optimization is applied to the cellular meta-material domain, using the architecture described in Section 5. The specifics of the training procedure and our use of the imitation learning technique DAGGER are explained in Section 6. Section 7 describes the software and hardware used. In Section 8, empirical evaluation demonstrates that composable energy surrogate models are able to solve cellular solid PDEs accurately and efficiently. Limitations and future work are discussed in Section 9.

Learning to optimize in collapsed bases
Solving PDEs like those that govern meta-material behavior can be framed as an optimization problem of finding a solution u which minimizes an energy E(u) subject to constraints. For mechanical meta-materials, E(u) is the stored elastic potential energy in the material. We propose a method for amortizing high-dimensional optimization problems where the objective has special conditional independence structure, such as that found in solving these PDEs. Consider the general problem of solving u * = arg min E(u) . (1) Here, u may be a vector in R d or may belong to a richer space of functions. Often we are only interested in either a small subset of the vector u * or the values the function u * takes on a small subdomain. To reflect this structure, we take the solution space to be the Cartesian product of a space of primary interest and a "nuisance" space. We denote the solutions as concatenations u = [x, y] where y is the object of interest, and x is the object whose value is not of interest to an application. x is roughly equivalent to auxiliary variables that often appear in probabilistic models, but that are marginalized away or discarded from the simulation. We can use this decomposition to frame Eq. 1 as a bi-level optimization problem: Consider the collapsed objective,Ẽ(y) = min x E(x, y). If it is possible to queryẼ(y) and its derivatives without ever representing x, we may perform collapsed optimization in the reduced basis of y, avoiding either performing optimization in the larger basis of u (Eq. 1), or performing bi-level optimization (Eq. 2) with an inner loop.
However,Ẽ is not usually available in closed form. We might consider approximatingẼ(y) via supervised learning. In the general case, this would require solv-ingẼ = min x E(x, y) for each example y we wish to include in our training set. This is the procedure used by many surrogate-based optimization techniques (Queipo et al., 2005;Forrester & Keane, 2009;Shahriari et al., 2015). The high cost of gathering each training example makes this prohibitive when x is high dimensional (such that the minimization is difficult) or when y is high dimensional (such that fitting a surrogate requires many examples).
In some cases, compositional structure in E may assist us with efficiently approximatingẼ. Many objectives may be represented as a sum: Given this decomposition, x i and x j are conditionally independent given y; i.e., if we constrain x i and y to take some values and perform minimization, the resulting x j or E j (x j , y) do not vary with the value chosen for x i . This follows from the partial derivative structure ∂Ei ∂xj = 0 for i = j. We propose to learn a collapsed objectiveẼ, which exploits conditional independence structure by represent-ingẼ(y) = iẼ i (y). This representation as a sum allows us to use min xi E i (x i , y) as targets for supervision, which may be found more cheaply than performing a full minimization. The learned approximations toẼ i may be composed to form an energy function with larger domain.
The language we use to describe this decomposition is intentionally chosen to reflect the conceptual similarity of our framework to collapsed variational inference (Teh et al., 2007) and collapsed Gibbs sampling (Geman & Geman, 1984;Liu, 1994). In these procedures, conditional independence is exploited to allow optimization or sampling to proceed in a collapsed space where nuisance random variables are marginalized out of the relevant densities. We exploit similar structure to these techniques, albeit in a deterministic setting. Other approaches to accelerating Eq. 2 which do not exploit (3) or directly modelẼ(y) include amortizing the inner optimization by predicting x * (y) = arg min x E(x, y) Figure 1. Cellular meta-materials. Top: materials at rest. Bottom: materials under axial compression, exhibiting periodic instability that varies with pore shape. The left two structures exhibit a negative Poisson's ratio, which does not occur in nature. (Kingma & Welling, 2013;Brock et al., 2017), or truncation of the inner loop, either deterministic (Wu et al., 2018;Shaban et al., 2018) or randomized to reduce bias (Tallec & Ollivier, 2017;Beatson & Adams, 2019).
The specific optimization procedure we accelerate is the numerical simulation of mechanical materials, where the objective corresponds to a physically meaningful energy, and the conditional independence structure arises from a spatial decomposition of the domain and the spatial locality of the energy density. We believe this spatial decomposition of the domain and energy could be generalized to learn collapsed energies for solving many other PDEs in reduced bases. This collapsed-basis approach may also be applicable to other bi-level optimization problems where the objective decomposes as a sum of local terms.

Mechanical meta-materials
Meta-materials are engineered materials with microstructure which results in macroscopic behavior not found in nature. They are often discussed in the context of materials achieving specific electromagnetic phenomena, such as negative refraction index solids and "invisibility cloaks" which electromagnetically conceal an object through engineered distortion of electromagnetic fields (Poddubny et al., 2013;Cai & Shalaev, 2010). However, they also hold great promise in other domains: mechanical meta-materials use substructure to achieve unusual macroscopic behavior such as negative Poisson's ratio and nonlinear elastic responses; pores and lattices undergo reversible collapse under large deformation, enabling the engineering of complex physical affordances in soft robotics (Bertoldi et al., 2017).
Meta-materials hold promise for modern engineering design, but are challenging to simulate as the microstructure necessitates a very fine finite element mesh, and as the nonlinear response makes them difficult to approximate with a macroscopic material model. Most work on meta-materials has thus relied on engineers and scientists to hand-design materials, rather than numerically optimizing meta-material substructure to maximize some objective (Ion et al., 2016).
We focus on building surrogate models for the twodimensional cellular solids investigated in Overvelde & Bertoldi (2014). These meta-materials consist of square "cells" of elastomer, each of which has a pore in its center. The pore shapes are defined by parameters α and β which characterize the pore shape in polar coordinates: r 0 is chosen such that the pore covers half the cell's volume: Constraints are placed on α and β to enforce a minimum material thicknesses and ensure that min θ r(θ) > 0 as in Overvelde & Bertoldi (2014).
These pore shapes give rise to complicated nonlinear elastic behavior, including negative Poisson's ratio and double energy wells (i.e., stored elastic energy which does not increase monotonically with strain). Realizations of this class of materials are shown under axial strain in Figure 1.
The continuum mechanics behavior of these elastomer metamaterials can be captured by a neo-Hookean energy model (Ogden, 1997). Let X ∈ R d , where d ≤ 3 in physical problems, be a point in the resting undeformed material reference configuration, and u(X) be the displacement of this point from reference configuration. The stored energy in a neo-Hookean solid is where W (u) is a scalar energy density over Ω. It is defined for material bulk and shear moduli µ and κ as: where F is the deformation gradient: The pore shapes influence the structure of these equations by changing the material domain Ω. The behaviour of these solids can be simulated by solving: where S = ∂W ∂F is known as the first Piola-Kirchoff stress, and where Eq. 9 defines a boundary condition. For example, G(u) = u − u b corresponds to a Dirichlet boundary condition; in our case, an externally imposed displacement. G(u) = ∂W ∂u − f b corresponds to an external force exerting a pressure on the boundary. Blue points locate control points of the splines used to represent the reduced-basis solutionũ.
To simulate these meta-materials, the PDE in Eq. 8 is typically solved via finite element analysis. Solving the PDEs arising from large mechanical meta-material structures is computationally challenging due to fine mesh needed to capture pore geometry and due to the highly nonlinear response induced by buckling under large displacements.
Solving the PDE in Eq. 8 corresponds to finding the u which minimizes the stored energy in the material subject to boundary conditions. That is, Eqs. 8 and 9 may be equivalently be expressed in an energy minimization form: In the next section, we use this energy-based perspective to facilitate learning energy surrogates which allow solving the PDE in a reduced basis of meta-material cell boundaries.

Composable energy surrogates
We apply the idea of learning collapsed objectives to the problem of simulating two-dimensional cellular mechanical meta-material behavior. The material response is determined by the displacement field u which minimizes the energy Ω W dX, subject to boundary conditions. We divide Ω into regular square subregions Ω i , which we choose to be cells with 2 × 2 arrays of pores, and denote the intersection of the subregion boundaries with B = ∂Ω 1 ∪ ∂Ω 2 ∪ . . . We let u i be the restriction of u to Ω i . We take the quantity of interest to be u B , the restriction of u to B, and the nuisance variables to be the restriction of u to Ω\B. The partitioning of Ω is shown in Figure 2.
The total energy in the material decomposes as a sum over the energy within each region: We introduce the collapsed energy of a component: This quantity is the lowest energy achievable by displacements of the interior of the cell Ω i , given the boundary conditions specified byũ i on ∂Ω i .Ẽ i (ũ i ) depends on the shape of the region Ω i , i.e., on the geometry of the pores. Rather than each possible pore shape having a unique collapsed energy function, we introduce the pore shape parameter ξ = (α, β) as an argument, replacingẼ i (ũ i ) withẼ(ũ i , ξ i ). The macroscopic behavior of the material is fully determined by this single collapsed energy func-tionẼ(ũ i , ξ i ). Given the true collapsed energy functions, we could accurately simulate material behavior in the reduced basis of the boundaries between each component Ω i . 1 We learn to approximate this collapsed energy function from data. This function may be duplicated and composed to simulate the material in the reduced basis B, an approach we term composable energy surrogates (CESs). A single CES is trained to approximate the functionẼ by fitting to su- where ξ i andũ i may be drawn from any distribution corresponding to anticipated pore shapes and displacements, and the targetsẼ(ũ i , ξ i ) are generated by solving the PDE in a small region Ω i with geometry defined by ξ i and withũ i imposed as a boundary condition. The CES may be duplicated and composed to approximate the total energy of larger cellular meta-materials.
Our ultimate goal is to efficiently solve for a reduced-basis displacement u B on B. Reduced-basis solving via CES may be cast as a highly-structured imitation learning problem. Consider using a gradient-based method to minimize the composed surrogate energy: whereÊ(ũ i , ξ i ) is the model's prediction ofẼ(ũ i , ξ i ), the collapsed energy of a single component. A sufficient condition for finding the correct minimum is for the "action" taken by the surrogate-the derivative of the energy approximation ∇ u BÊ -to match the "action" taken by an expert-the total derivative, ∇ u B min u / ∈B E(u)-along the optimization trajectory. If so, the surrogate will follow the trajectory of a valid, if non-standard, bilevel gradient-based procedure for minimizing the energy, corresponding to (2). Given an imperfect surrogate, the error in the final solution could trivially be bounded in terms of the error in approximating ∇ u B min u / ∈B E(u) by ∇ u BÊ along the trajectory, and the number of gradient steps taken. This observation informs our model, training, and data collection procedures, described in the following sections.

Model architecture
Our CESs take the form of a neural architecture, designed to respect known properties of the true potential energy and to maximize usefulness as surrogate energy to be minimized via a gradient-based procedure. The effects of these design choices are quantified via an ablation study in the appendix.
Reduced-basis parameterization. We require a vector representation for the functionũ. We use one cubic spline for each horizontal and vertical displacement function along each face of the square, with evenly spaced control points and "not-a-knot" boundary conditions. Our vector representation is u ∈ R 2n , formed from the horizontal and the vertical displacement values at each of the n control points. Splines on adjacent faces share a control point at the corner. Using N control points to parameterize the function along each face requires n = 4 * (N − 1) control points to parameterize a 1d function around a single cell. For all experiments we use N = 10 control points along each edge, resulting in a vector u with 72 entries.
Model structure and loss. Our surrogate energy model is: where f φ is a neural network with parameters φ and R removes rigid-body rotation and translation. Our loss function is L = L 0 + L 1 + L 2 , which is the sum of losses on the 0th, 1st and 2nd energy derivatives: Log-stiffness loss Cosine distance between Hessian-vector products v ∼ N (0, I 2n ) Projection vector for Hessian .
In the above, ∇ u and ∇ 2 u are the gradient and Hessian of the surrogate energyÊ or the ground-truth energyẼ with respect to u, and v is sampled independently for each training example in a batch.
Invariance to rigid body transforms. The true stored elastic energy is invariant to rigid body transformation of a solid. This invariance may be hard for a neural network to learn exactly from data. We define a module R which applies Procrustes analysis, a procedure that finds and applies the rigid body transformation which minimizes the Euclidean distance to a reference, for which we use the rest configuration. This is differentiable and closed-form.
Encoding a linear elastic bias. The energy is approximated well by a linear elastic model when at rest:Ẽ i (ũ i ) ≈ R(u i ) T A i R(u i ) for a stiffness matrix A i depending on ξ i . We scale our net's outputs by ||R(u i )|| 2 2 so that it needs only capture a "scalar stiffness" E /||R(ui)|| 2 2 accounting for the geometry of A i given ξ i and for deviation from the linear elastic model.
Parameterizing the log-stiffness. The energy of a compo-nentẼ i (u 0,i ) is nonnegative, and the ratio of energy to a linear elastic approximation varies over many orders of magnitude. We thus parameterize the log of the scalar stiffness with our neural network f φ rather than the stiffness.
Log-stiffness loss. We wish to find neural network parameters φ which lead to accurate energy predictions for many different orders of magnitude of energy and displacement. Minimizing the 2 loss between predicted and true energies penalizes errors in predicting large energies more than proportional errors predicting small energies. Instead, we take the 2 loss between the predicted logstiffness f φ (R(u), ξ) and the effective ground-truth logstiffness, logẼ (ũ) /||R(u)|| 2 2 . Sobolev training with gradients. When derivatives of a target function are available, training a model to match these derivatives ("Sobolev training") can aid generalization (Czarnecki et al., 2017). Accuracy of CES' gradients is crucial to an accurate solution trajectory. We obtain groundtruth gradients cheaply via the adjoint method (Lions, 1971). Given a solution u i to the PDE in Ω i with boundary conditionsũ i , the gradient ∇ũ iẼi (ũ i ) requires one linear system solve, with the same cost as one Newton step while solving the PDE (Mitusch et al., 2019). The spline is a linear map M from u i toũ i in the finite element basis. We can thus efficiently compute ∇ uiẼi (ũ i ) = M T ∇ũ iẼi (ũ i ). The gradient of our model, ∇ uiÊφ (u i , ξ i ), may be computed with one backward pass.
Sobolev training with Hessian-vector products. Given a solution and gradient, computing ∇ 2 uẼ requires 2n linear system solves-one for each entry in u. As 2n = 72 is much smaller than the number of Newton steps we need to solve the PDE for moderate displacements, we expect the increased fidelity from a 2nd-order approximation of the energy to be worth this added computation. Computing the full Hessian of the surrogate energy, ∇ 2 uiÊφ (u i , ξ i ), would require 2n backward passes. Instead we train on Hessianvector products, which require only a single backward pass additional to that required for the gradient.
Cosine distance loss for Sobolev training. Ground-truth gradient and Hessian values vary over many orders of magnitude, roughly corresponding to lower and higher energy displacements. We wish our model to be robust to outliers and accurate across a range of different conditions. Rather than placing an 2 loss on the gradient and Hessian-vector products as in Czarnecki et al. (2017), we minimize the cosine distance between ground truth and approximate gradients and Hessians, which is naturally bounded in [0, 1].

Data and training
We carry out data collection in two phases. First, we collect training and validation datasets using Hamiltonian Monte Carlo (Duane et al., 1987) to preferentially sample displacements which correspond to lower energy modes. Next, we perform dataset aggregation (Ross et al., 2011) to augment the dataset such that the learned energy model will be accurate on the states it encounters when deployed.

Solving the PDE
To collect training data, we use the reduced-basis displacementũ corresponding to a vector of spline coefficients u as the boundary condition around a domain Ω representing a 2 × 2-pore subdomain, and solve the PDE Div F (u) = 0 ∈ Ω u =ũ ∈ ∂Ω using a load-stepped relaxed Newton's method (Sheng et al., 2002). The relaxed Newton's method takes the iteration: Above, 0 < λ < 1 is the relaxation parameter (analogous to a step size), and u is the vector of coefficients defining u in the FEA basis. Newton's method requires an initial guess which is sufficiently close to the true solution (Kythe et al., 2004). Smaller relaxation parameters yield a greater radius of convergence but necessitate more steps to solve the PDE.
The radius of convergence can also be aided by loadstepping: solving the PDE for a sequence of boundary conditions, annealing from an initial boundary condition for which we have a good initial guess (e.g., the rest configuration) to a final boundary conditionũ, using the solution to the previous problem as an initial guess for Newton's method for the next problem. We find that combining load stepping with a relaxed Newton's method allows problems to be solved more efficiently than using either alone. Except where specified, we linearly anneal from rest toũ over 10 load steps and use a relaxation parameter λ = 0.1.

Initial dataset collection
We wish to train on a wide variety of displacement boundary conditions. Solution procedures minimize the energy: thus, lower energy modes will be encountered in the solve. We choose a shaping distribution where the density is the product of a Boltzmann density exp{Ẽ}/Z and a Gaussian density N (x(u);μ, Σ), wherex(u) ∈ R 2×2 is a macroscopic strain tensor 2 represented by u,μ is a target macroscopic strain drawn from an i.i.d. Gaussian with standard deviation 0.15, and Σ is set to (μ •μ) −1 .
Given a solution to the PDE, the log-density and its displacement may be cheaply computed (the latter via the adjoint method). Making use of these gradients, we sample data points with Hamiltonian Monte Carlo (HMC). After sampling a data point, we compute the corresponding Hessian and save the tuple (u, ξ,Ẽ, ∇ uẼ , ∇ 2 uẼ ) as a data point. We initialize each HMC data collector by sampling a macroscopic displacement target and a random pore shape. We do not use load-stepping, instead using the solution for the u used in the previous iteration of HMC's leapfrog integration as an initial guess for solving the PDE. We randomize HMC hyperparameters for each collector to attempt to minimize the impact of specific settings: see the appendix for exact ranges. We sample 55,000 training examples and 5000 validation examples altogether. We visualize displacements drawn from this distribution in the appendix.

Data aggregation
The procedure of deploying the surrogate defies standard i.i.d. assumptions in supervised machine learning. That is, when deployed, the surrogate encounters states determined by the energy it defines and on the boundary conditions placed on the composed body. Depending on this energy and on the boundary conditions, the surrogate may encounter states which do not resemble those sampled with HMC. This problem-that training an agent to predict expert actions with supervised learning leads to trajectories dissimilar to those on which it was trained-is a central concern in the imitation learning literature. A number of solutions exist (Schroecker & Isbell, 2017). One such technique is dataset aggregation, or DAGGER (Ross et al., 2011), an extension of earlier approaches SEARN (Daumé et al., 2009) and SMILe (Ross & Bagnell, 2010), which reduces imitation learning or structured prediction to online learning.
In DAGGER, a policy is deployed and trajectories are collected. The expert is queried for ground-truth actions on the states in these trajectories. The state-action pairs are appended to the dataset, and the policy is retrained on this dataset. This process of deployment, querying, appending data, and retraining, is iterated. Under appropriate assumptions, the instantaneous regret of the learned policy vanishes with the number of iterations, i.e., the learned policy will match the expert policy on its own trajectories. Ross et al. (2011) present DAGGER as a method for discrete action spaces. Our problem has a continuous action space: the gradient of the energy in a cell. We do not investigate whether it is possible to generalize DAGGER's regret guaran- tees to continuous action spaces, but the intuition holds that we wish our model to "imitate" the finite element "expert" on the optimization trajectories the model produces.
We initialize our training data as described in the preceding section. We then apply DAGGER by iterating: (i) training the surrogate; (ii) composing surrogates and finding the displacements which minimize the composed energy; (iii) sampling displacements which lie along the surrogate's solution path, querying the ground-truth energy and energy derivatives using the finite element model, and adding these new data points to the dataset. We visualize displacements generated by DAGGER in the appendix.
We use Ray (Moritz et al., 2018) to run distributed workloads on Amazon EC2. The initial dataset is collected using 80 M4.xlarge CPU spot workers. While training the surrogate, we use a GPU P3.large driver node to train the model, and 80 M4.xlarge CPU spot worker nodes performing DAG-GER in parallel. These workers receive updated surrogate model parameters, compose and deploy the surrogate, sample displacements along the solution path, query the finite element model for energy and derivatives, and return data to the driver node. Initial dataset collection and model training with DAGGER each take about one day in wall-clock time.

Empirical evaluation
In this section we demonstrate the ability of Composable Energy Surrogates (CES) to efficiently produce accurate solutions. We compare wall-clock computation time and solution accuracy of CES to that of FEA at varying fidelities.
We consider the systems constructed in Overvelde & Bertoldi (2014): structures with an 8 × 8 array of pores, corresponding to a 4 × 4 assembly of our surrogates, each of which represents a 2 × 2-pore component. We sample pore shapes from a uniform distribution over the valid shapes defined in Overvelde & Bertoldi (2014). For DAGGER, we sample macroscopic vertical axial strain magnitudes from U(0., 0.3), and choose to apply compression with probability 0.8 (as compressive displacements involve more interesting pore collapse) or tension with probability 0.2.
We compare our composed surrogates to finite element analysis with a range of different-fidelity meshes under axial compression and tension with a macroscopic displacement of 0.125L 0 , where L 0 is the original length of the solid. See the appendix for details of the finite element meshes. Comparison is carried out for seven pore shapes: ξ = (0, 0), corresponding to circular pores, and six ξ sampled from a uniform distribution over pore parameters defined as valid in Overvelde & Bertoldi (2014) via rejection sampling.
For the composed surrogate, we use PyTorch's L-BFGS routine to minimize the energy, with fixed step size 0.25 and PyTorch's default criteria for checking convergence. We attempt to solve each finite element model with Newton's method with load steps [1,2,5,10,20] and relaxation parameters [0.9, 0.7, 0.4, 0.1, 0.05]. We record the time taken for the fastest solve which converges. Under compression these solids exhibit nonlinear behavior, and only the more conservative solves converge. Under tension they behave closer to a linear elastic model, and Newton's method converges quickly. Measurements are taken on an Amazon AWS M4.xlarge EC2 CPU instance. Using a GPU could provide further acceleration for the composed surrogate.
We measure error in the solution and in the macroscopic energy. The former is ||û − u * || 2 2 , whereû and u * are the approximation and ground-truth evaluated at the spline control points. We also measure the relative error, |Ê(û)−E * (u * )| /E * (u * ), whereÊ(û) is the approximated energy of the approximate solution, and E * (u * ) is the groundtruth energy of the ground-truth solution. As the energy function determines macroscopic behavior, accuracy of this energy is a potential indicator of a model's ability to generalize to larger structures. The highest-fidelity finite element model is used for the ground truth E * and u * , and has an error of zero on both metrics by definition. Multiple solutions exist as energy is preserved under vertical and horizontal flips, so before comparing a solutionû to the ground-truth u * we programatically check these flips and use the flip which minimizes the Euclidean distance. Figure 3 shows our evaluation. Composed energy surrogates are more efficient than high-fidelity FEA simulations yet more accurate than low-fidelity simulations. CES produces solutions with equivalent 2 error to FEA solutions which need an order of magnitude more variables and computation time, or produces solutions with an order of magnitude less 2 error than FEM solutions requiring the same computation. This gap increases to several orders of magnitude when we consider percentage error in the predicted strain energy. We visualize the ground-truth and the CES approximation in Figure 4. See the appendix for visualization of the FEM solutions and of CES for the remaining structures.

Limitations and opportunities
Use of DAGGER. We use DAGGER to assist CES to match the ground-truth on the states encountered during the solution trajectory. DAGGER requires one to specify in advance the conditions under which the surrogate will be deployed, limiting the surrogate's use for arbitrary downstream tasks. Investigating CES's ability to generalize to novel deployment conditions-and designing surrogates which can do so effectively-is an important direction for future work.
Error estimation, refinement, and guarantees. Finite element methods enjoy advantageous properties. As one uses a higher-resolution mesh or higher-order elements, the solution in the finite element basis approaches the true solution. This provides a straightforward way to estimate the error (compare to the solution in a more-refined basis) and control it (via refinement). CES currently lacks these properties.
Finite element baseline. There is an immense body of work on finite element methods and iterative solvers. We try to provide a representative baseline, but our work should not be taken as a comparison with the "state-of-the-art". We aim to show that composable machine-learned energy surrogates enjoy some advantages over a reasonable baseline, and hold promise for scalable amortization of solving modular PDEs.
Hyperparameters. Both our method and the finite element baseline rely on a multitude of hyperparameters: the size of the spline reduced basis; the size and learning rate of the neural network; the size and degree of the finite element approximation; and the specific variant of Newton's method to solve the finite element model. We do not attempt a formal, exhaustive search over these parameters.
Known structure. We leave much fruit on the vine in terms of engineering structure known from the true equations into our surrogate. For example, we do not engineer invariance of the energy to flips and rotations of the spline coefficients. One could also use a more expressive normalizer than ||u|| 2 2 , e.g. the energy predicted by a coarse-grained linear elastic model, or exploit spatially local correlation, e.g. by using a 1-d convolutional network around the boundary of the cell.

Conclusion
We present a framework for collapsing optimization problems with a local bilevel structure by learning composable energy surrogates. This framework is applied to amortizing the solution of PDEs corresponding to mechanical metamaterial behavior. Learned composable energy surrogates are more efficient than high-fidelity FEA yet more accurate than low-fidelity FEA, occupying a new point on the Pareto frontier. We believe that learning composable energy surrogates could accelerate metamaterial design, as well as design and identification of other systems described by PDEs with parametric modular structure.

Acknowledgements
We would like to thank Alexander Niewiarowski for numerous helpful discussions about continuum mechanics and FEA, Ari Seff for help finding a particularly difficult bug in our data pipeline, and Maurizio Chiaramonte for inspiring early conversations about metamaterials and model order reduction. This work was funded by NSF IIS-1421780 and by the Princeton Catalysis Initiative.
bounding circle divided by the resolution. It should be noted that cells can be much smaller than this maximum size, or there can be many more cells than the resolution parameter would imply, as mshr will place one cell vertex on each point used to construct the domain geometry (i.e. each point in the pore polygon).

Benchmark visualizations
In the following pages we visualize the solutions found for each pore by each FEA mesh and by CES. For each pore we use six different finite element meshes. These respectively used [4,8,16,32,48,64] points used to define the geometry of each pore, and minimum of [1,2,4,8,12,16] internal mesh vertices along a given axis per pore. Given these parameters and the geometry of the material domain, meshes were created using the automatic mesh generation tool from mshr (the mesh generation component of FEniCS). We include the number of degrees of freedom in the finite element basis in each plot. We superimpose the solution found with CES in red dots on the solution found with FEA. The CES solution has 690 degrees of freedom in all cases.