<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Using scipy iterative solvers
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.
Step1:  Generating the microstructure
The microstructure is generated by means of the gen_ellipsoid.py script, which is listed below.
Step2:  The above script should be invoked as follows
python gen_ellipsoid.py tutorial.json
where the tutorial.json file holds the geometrical parameters.
Step3:  The resulting microstructure is to be saved into the example.h5 file.
Step4:  The microstructure is then retrieved as follows.
Step5:  We can retrieve dim and sym from the dimensions of the microstructure.
Step6:  And we can check visually that everything went all right.
Step7:  Creating the basic objects for the simulations
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.
Step8:  We then define instances of the IsotropicLinearElasticMaterial class for all three materials.
Step9:  We want to solve the Lippmann–Schwinger equation, which reads
\begin{equation}
\bigl(\mathbf C-\mathbf C_{\mathrm{ref}}\bigr)^{-1}
Step10:  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}
Step11:  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!).
Step12:  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.
Step13:  We know that the linear operator $A$ is definite. We can therefore use the conjugate gradient method to solve $A\cdot x=b$.
Step14:  The resulting solution, $x$, must be reshaped into a $N_1\times\cdots\times N_d\times s$ array.
Step15:  We can plot the $\tau_{xy}$ component.
Step16:  And compute the associated strain field.
Step17:  And plot the $\varepsilon_{xy}$ component.
<ASSISTANT_TASK:>
Python Code:
import h5py as h5
import matplotlib.pyplot as plt
import numpy as np
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
from scipy.sparse.linalg import cg, LinearOperator
%matplotlib inline
plt.rcParams['figure.figsize'] = (12, 8)
Explanation: Using scipy iterative solvers
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.
End of explanation
%pfile gen_ellipsoid.py
Explanation: Generating the microstructure
The microstructure is generated by means of the gen_ellipsoid.py script, which is listed below.
End of explanation
%pfile tutorial.json
Explanation: The above script should be invoked as follows
python gen_ellipsoid.py tutorial.json
where the tutorial.json file holds the geometrical parameters.
End of explanation
%run gen_ellipsoid.py tutorial.json
Explanation: The resulting microstructure is to be saved into the example.h5 file.
End of explanation
with h5.File('./tutorial.h5', 'r') as f:
    phase = np.asarray(f['phase'])
    
Explanation: The microstructure is then retrieved as follows.
End of explanation
dim = phase.ndim
sym = (dim*(dim+1))//2
Explanation: We can retrieve dim and sym from the dimensions of the microstructure.
End of explanation
plt.imshow(phase);
Explanation: And we can check visually that everything went all right.
End of explanation
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
Explanation: Creating the basic objects for the simulations
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.
End of explanation
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)
Explanation: We then define instances of the IsotropicLinearElasticMaterial class for all three materials.
End of explanation
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)
Explanation: We want to solve the Lippmann–Schwinger equation, which reads
\begin{equation}
\bigl(\mathbf C-\mathbf C_{\mathrm{ref}}\bigr)^{-1}:\boldsymbol\tau+\boldsymbol\Gamma_{\mathrm{ref}}[\boldsymbol\tau]=\mathbf E,
\end{equation}
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
\begin{equation}
\bigl(\mathbf C^h-\mathbf C_{\mathrm{ref}}\bigr)^{-1}:\boldsymbol\tau^h+\boldsymbol\Gamma_{\mathrm{ref}}^h[\boldsymbol\tau^h]=\mathbf E,
\end{equation}
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
\begin{equation}
A\cdot x=b,
\end{equation}
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).
End of explanation
a = MyLinearOperator(phase, C_m, C_i, C_ref)
Explanation: 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.
Solving the Lippmann–Schwinger equation
We are now ready to solve the equation. We first create an instance of the linear operator $A$.
End of explanation
eps_macro = np.zeros((sym,), dtype=np.float64)
eps_macro[-1] = np.sqrt(2.)
Explanation: 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!).
End of explanation
b_arr = a.empty_arr()
b_arr[...] = eps_macro
b = b_arr.ravel()
Explanation: 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.
End of explanation
x, info = cg(a, b)
assert info == 0
Explanation: We know that the linear operator $A$ is definite. We can therefore use the conjugate gradient method to solve $A\cdot x=b$.
End of explanation
tau = x.reshape(a.arr_shape)
Explanation: The resulting solution, $x$, must be reshaped into a $N_1\times\cdots\times N_d\times s$ array.
End of explanation
plt.imshow(tau[..., -1]);
Explanation: We can plot the $\tau_{xy}$ component.
End of explanation
eps = a.tau2eps.apply(tau)
Explanation: And compute the associated strain field.
End of explanation
plt.imshow(eps[..., -1]);
plt.plot(eps[63, :, -1])
Explanation: And plot the $\varepsilon_{xy}$ component.
End of explanation
<END_TASK>