The aim of this notebook is to show how the scipy.sparse.linalg module can be used to solve iteratively the Lippmann–Schwinger equation.
The problem at hand is a single ellipsoid in a periodic unit-cell, subjected to a macroscopic strain $\mathbf E$. Again, we strive to make the implementation dimension independent. We therefore introduce the dimension dim of the physical space, and the dimension sym = (dim*(dim+1))//2 of the space of second-order, symmetric tensors.
We start by importing a few modules, including h5py, since input and output data are stored in HDF5 format.
In [ ]:
import h5py as h5
import matplotlib.pyplot as plt
import numpy as np
In [ ]:
import janus
import janus.material.elastic.linear.isotropic as material
import janus.operators as operators
import janus.fft.serial as fft
import janus.green as green
In [ ]:
from scipy.sparse.linalg import cg, LinearOperator
In [ ]:
%matplotlib inline
In [ ]:
plt.rcParams['figure.figsize'] = (12, 8)
The microstructure is generated by means of the gen_ellipsoid.py script, which is listed below.
In [ ]:
%pfile gen_ellipsoid.py
The above script should be invoked as follows
python gen_ellipsoid.py tutorial.jsonwhere the tutorial.json file holds the geometrical parameters.
In [ ]:
%pfile tutorial.json
The resulting microstructure is to be saved into the example.h5 file.
In [ ]:
%run gen_ellipsoid.py tutorial.json
The microstructure is then retrieved as follows.
In [ ]:
with h5.File('./tutorial.h5', 'r') as f:
phase = np.asarray(f['phase'])
We can retrieve dim and sym from the dimensions of the microstructure.
In [ ]:
dim = phase.ndim
sym = (dim*(dim+1))//2
And we can check visually that everything went all right.
In [ ]:
plt.imshow(phase);
We first select the elastic properties of the inclusion, the matrix, and the reference material. For the latter, we select a material which is close to the matrix, but not equal, owing to the $(\mathbf C_{\mathrm m}-\mathbf C_{\mathrm{ref}})^{-1}$ factor in the Lippmann–Schwinger equation.
In [ ]:
mu_i, nu_i = 10., 0.3 # Elastic properties of the ellipsoidal inclusion
mu_m, nu_m = 1., 0.2 # Elastic properties of the matrix
mu_ref, nu_ref = 0.99*mu_m, nu_m # Elastic properties of the reference material
We then define instances of the IsotropicLinearElasticMaterial class for all three materials.
In [ ]:
C_i = material.create(mu_i, nu_i, dim=dim)
C_m = material.create(mu_m, nu_m, dim=dim)
C_ref = material.create(mu_ref, nu_ref, dim=dim)
type(C_i)
We want to solve the Lippmann–Schwinger equation, which reads
where $\mathbf C=\mathbf C_{\mathrm i}$ in the inclusion, $\mathbf C=\mathbf C_{\mathrm m}$ in the matrix, and $\boldsymbol\Gamma_{\mathrm{ref}}$ is the fourth-order Green operator for strains. After suitable discretization, the above problem reads
where $\mathbf C^h$ denotes the local stiffness, discretized over a cartesian grid of size $N_1\times\cdots\times N_d$; in other words, it can be viewed as an array of size $N_1\times\cdot\times N_d\times s\times s$ and $\boldsymbol\Gamma_{\mathrm{ref}}^h$ is the discrete Green operator. The unknown discrete polarization field $\boldsymbol\tau^h$ ($N_1\times\cdots\times N_d\times s$ array) is constant over each cell of the cartesian grid. It can be assembled into a column vector, $x$. Likewise, $\mathbf E$ should be understood as a macroscopic strain field which is equal to $\mathbf E$ in each cell of the grid; it can be assembled into a column vector, $b$.
Finally, the operator $\boldsymbol\tau^h\mapsto\bigl(\mathbf C^h-\mathbf C_{\mathrm{ref}}\bigr)^{-1}:\boldsymbol\tau^h+\boldsymbol\Gamma_{\mathrm{ref}}^h[\boldsymbol\tau^h]$ is linear in $\boldsymbol\tau^h$ (or, equivalently, $x$); it can be assembled as a matrix, $A$. Then, the discrete Lippmann–Schwinger equation reads
which can be solved by means of any linear solver. However, two observations should be made. First, the matrix $A$ is full; its assembly and storage might be extremely costly. Second, the matrix-vector product $x\mapsto A\cdot x$ can efficiently be implemented. This is the raison d'être of a library like Janus!
These observation suggest to implement $A$ as a linearOperator, in the sense of the scipy library (see reference).
In [ ]:
class MyLinearOperator(LinearOperator):
def __init__(self, phase, C_m, C_i, C_ref):
dim = phase.ndim
sym = (dim*(dim+1))//2
alpha_i = 1./dim/(C_i.k-C_ref.k)
beta_i = 1./2./(C_i.g-C_ref.g)
alpha_m = 1./dim/(C_m.k-C_ref.k)
beta_m = 1./2./(C_m.g-C_ref.g)
T = np.array([operators.isotropic_4(alpha_i, beta_i, dim),
operators.isotropic_4(alpha_m, beta_m, dim)])
self.tau2eps = operators.block_diagonal_operator(T[phase])
self.green = green.filtered(C_ref.green_operator(),
phase.shape, 1.,
fft.create_real(phase.shape))
self.arr_shape = phase.shape+(sym,)
n = np.product(self.arr_shape)
super().__init__(np.float64, (n, n))
def _matvec(self, x):
tau = x.reshape(self.arr_shape)
eta = np.zeros_like(tau)
self.tau2eps.apply(tau, eta)
eta += self.green.apply(tau)
y = eta.ravel()
return y
def _rmatvec(self, x):
return self._matvec(x)
def empty_arr(self):
return np.empty(self.arr_shape)
The constructor of the above class first computes the local map $\boldsymbol\tau^h\mapsto\bigl(\mathbf C^h-\mathbf C_{\mathrm{ref}}\bigr)^{-1}$. Then, it implements the operator $\boldsymbol\tau^h\mapsto\bigl(\mathbf C^h-\mathbf C_{\mathrm{ref}}\bigr)^{-1}:\boldsymbol\tau^h$. The resulting operator is called tau2eps.
The constructor also implements a discrete Green operator, associated with the reference material. Several discretization options are offered in Janus. The filtered Green operator is a good option. TODO Use the Willot operator instead.
Finally, the operator $A$ is implemented in the _matvec method, where attention should be paid to the fact that x is a column-vector, while green and tau2eps both operates on fields that have the shape of a symmetric, second-order tensor field defined over the whole grid, hence the reshape operation. It is known that the operator $A$ is symmetric by construction. Therefore, the _rmatvec method calls _matvec.
We are now ready to solve the equation. We first create an instance of the linear operator $A$.
In [ ]:
a = MyLinearOperator(phase, C_m, C_i, C_ref)
We then create the macroscopic strain $\mathbf E$ that is imposed to the unit-cell. In the present case, we take $E_{xy}=1$ (beware the Mandel–Voigt notation!).
In [ ]:
eps_macro = np.zeros((sym,), dtype=np.float64)
eps_macro[-1] = np.sqrt(2.)
We then populate the right-hand side vector, $b$. b_arr is the column-vector $b$, viewed as a discrete, second order 2D tensor field. It is then flattened through the ravel method.
In [ ]:
b_arr = a.empty_arr()
b_arr[...] = eps_macro
b = b_arr.ravel()
We know that the linear operator $A$ is definite. We can therefore use the conjugate gradient method to solve $A\cdot x=b$.
In [ ]:
x, info = cg(a, b)
assert info == 0
The resulting solution, $x$, must be reshaped into a $N_1\times\cdots\times N_d\times s$ array.
In [ ]:
tau = x.reshape(a.arr_shape)
We can plot the $\tau_{xy}$ component.
In [ ]:
plt.imshow(tau[..., -1]);
And compute the associated strain field.
In [ ]:
eps = a.tau2eps.apply(tau)
And plot the $\varepsilon_{xy}$ component.
In [ ]:
plt.imshow(eps[..., -1]);
In [ ]:
plt.plot(eps[63, :, -1])