Alejandro Queiruga's Presentation

## Numerical Experiments of State-based  Peridynamics and Influence Functions in Classical Problems
### Alejandro F. Queiruga  
### George J. Moridis
#### June 7, 2018
#### Lawrence Berkeley National Lab

---

## Bottom Line Upfront

Peridynamics does not pass tests for a mechanics discretization:
1. It does not represent constant strain modes exactly.
2. The only converging formulation is unstable, especially around fracture tips.
3. We can fix the instability, but that doesn't make #1 exact.

`git clone` my implementation PeriFlakes:
[https://github.com/afqueiruga/PeriFlakes](https://github.com/afqueiruga/PeriFlakes)  
The Jupyter notebooks replicate the graphs in this presentation.

Today, I'm going to talk about the constant strain modes.

---

## Outline

- The problem I was trying to solve
- The simulation I managed to program
- How do we analyze this method?
- PeriFlakes
- Constant strain problems
- Boundary conditions
- Briefly: Linear fracture mechanics

---

## Motivating Problem

Shale natural gas, Enhanced geothermal, Escape of gasses from nuclear
waste,  Bubbling of methane through oceanic sediments

---

## Peridynamic Representation

---

## Computational Framework

---

## Simulation Results

---

## Natural Fracture Interaction

---

## Videos are not results

It looks good, but not under actual scrutiny:

>Queiruga, A. F. and G. J. Moridis, "Numerical experiments on 
 > the convergence properties of state-based peridynamic laws and 
>  influence functions in two-dimensional problems." Computer 
 > Methods in Applied Mechanics and Engineering 322 (2017): 
 > 97-122.

**It is critical to verify and validate numerical methods on known
solutions and expected discretization refinement behavior.**

- Can't jump to dynamics and fracturing simulations immediately.

---

## What is Peridynamics?

- We can describe linear mechanics as an operator $L$ in a continuum.
- PD defines an operator that depends on a damage $\alpha$:
\begin{equation}
\rho \ddot{\mathbf{u}} = L(\alpha)\mathbf{u} + \mathbf{b}
\end{equation}
with evolution laws for the damage,
\begin{equation}
\dot{\alpha} = \min{\left(f(\alpha,\mathbf{u}), 0\right)}
\end{equation}
- Only considering static 2D plane strain linear elasticity:
\begin{equation}
0 = L_{\alpha=1} \mathbf{u}
\end{equation}
- Considering it only as a discretized operator.

What the properties of the Peridynamics specification of $L(\alpha)$?

---

## How do we know our codes are right?

- Unit tests against certain solutions
- Convergence to known analytical solutions
- Convergence to manufactured solutions
- Convergence to super-refined self solutions

For a given method, some problems are in an exactly-representable
solution space. E.g., Quadratic finite elements solve quadratic solutions exactly.

Which solutions are exactly representable by Peridynamics?

---

## Constant Strain Problems

- Three basic modes of deformation with a constant gradient
- Fundamental solutions from basic mechanics
- FEM is correct to within machine precision
- This is the patch test.

These are **Unit Tests** in **all** of my codes! (PD, FEM, FDM, RKPM, etc.)
```python
self.AssertTrue( error < 1.0e-12 )
```

---

## Thin Pressurized Crack

Also check for convergence to more complicated analytical solutions.

<img src="assets/figures/crack_diagram.png" style="float:left; height: 60vh" >
\begin{equation}
\hat{Z}\left(z = x+iy\right)=P\left(\sqrt{z^{2}-L^{2}}-z\right)
\end{equation}

\begin{equation}
u\_{x} = 
\frac{1}{2G}\left(\frac{2-4\nu}{2}\mathrm{Re}\left[\hat{Z}\right]-y\mathrm{Im}\left[Z\right]\right)
\end{equation}
\begin{equation}
u\_{y} = \frac{1}{2G}\left(\frac{4-4\nu}{2}\mathrm{Im}\left[\hat{Z}\right]-y\mathrm{Re}\left[Z\right]\right)
\end{equation}

Calculate an error against the
analytical solution $u^a$ to judge the numerical solution $u^N$:
\begin{equation}
e = \frac{ \int \left\|u^N-u^a \right\|^2 \mathrm{d}x }{ \int
\left\|u^a \right\|^2 \mathrm{d}x }
\end{equation}

---

## Fixed-$\delta$ Convergence

Add more grid points keeping $w(r;\delta)$ fixed.

---

### The support size grows with $m$

<img src="assets/plots/peri_support_size.png" style="max-height: 50vh" >
*Periflakes/computation.ipynb*
Larger $\delta$ vs $h$ means more points interact.

---

### This means the matrix fills up...

<img src="assets/plots/sparsity_1.5.png" style="max-height: 50vh" >
\begin{equation}
m=1.5
\end{equation}
*Periflakes/computation.ipynb*

---

### This means the matrix fills up...

<img src="assets/plots/sparsity_2.5.png" style="max-height: 50vh" >
\begin{equation}
m=2.5
\end{equation}
*Periflakes/computation.ipynb*

---

### This means the matrix fills up...

<img src="assets/plots/sparsity_3.5.png" style="max-height: 50vh" >
\begin{equation}
m=3.5
\end{equation}
*Periflakes/computation.ipynb*

---

### This means the matrix fills up...

<img src="assets/plots/sparsity_4.5.png" style="max-height: 50vh" >
\begin{equation}
m=4.5
\end{equation}
Even if we don't assemble a matrix, this still shows data access patterns
of the program.

### The number of nonzeros grows:

<img src="assets/plots/peri_nnz_growth.png" style="max-height: 50vh" >
Even if we don't assemble a matrix, this still shows data access patterns
 of the program.
*Periflakes/computation.ipynb*

---

## Fixed-$m$ Convergence

- Changing $h$ but not $\delta$ changes the stencil.
- Refining $h$ but not $\delta$ leads to intractable computation!

- By keeping $m = \frac{\delta}{h}$ constant, the stencil stays the same.
- $m$ is the discretization-independent (hyper)parameter of the numerical method.

---

## Hyperparameters

Our goal is to **obtain a computer program** that solves the problems we
want accurately and efficiently.

In the Peridynamics scope, we have the following hyperparameters:

1. Formulation $\mathbf{T}\left[x\right]\left<\xi\right>$- discrete choice
2. Influence function $w\left(r\right)$ - discrete choice
3. Relative support radius, $m=\delta/h$ - continuous
4. Boundary condition - discrete choice
5. Stabilizer parameter $A$ - continuous
6. Smoothing function $w^s$ - discrete choice (including None)

These are **arbitrary**. We need to determine which combinations are
acceptable, and which of those are the best.

---

## Variety of Formulations

It is well known that bond-based does not work due to $\nu$ problem.

The are three classes of state-based models:

1. The ordinary dilation-state methods
  1. Silling's
  2. Oterkus's
2. The deformation gradient (F) correspondence methods
  1. Silling's original formulation
  1. Littlewood 2010/11's stabilization
  2. Silling 2017's stabilization (back-to-back in CMAME to Queiruga, 2017!)
  4. Post-processing smoothing
3. Non-ordinary methods
  - I haven't tried any yet.

## Smoothing

My ad hoc approach was to post process the solution:
\begin{equation}
\mathbf{u}^s = \frac{\int\_\mathcal{H} w^s \alpha \mathbf{u} \mathrm{d}x'}{\int\_\mathcal{H} w^s \alpha \mathrm{d}x'}
\end{equation}

The smoothing kernel $w^s$ doesn't have to be the same $w$.

This worked to stabilize the fracture path in my original hydraulic fracturing simulations,
but the method never passed any other validations!

## The Stabilizers

We can add a stabilizer to the force:
\begin{equation}
\mathbf{T}\left<\xi\right> = w\alpha \sigma \mathbf{K}^{-1} \mathbf{\xi} + \mathbf{T}^s
\end{equation}

Silling's 2017 paper:
\begin{equation}
\mathbf{T}^s = w\alpha A \frac{18 K}{ w_0 \pi \delta^5} \left( \mathbf{Y} - \mathbf{F}\xi \right),
\quad w\_0 = \int_\mathcal{H} w\alpha \mathrm{d}x'
\end{equation}

Littlewood's 2010 paper:
\begin{equation}
\mathbf{T}^s = w\alpha A \frac{18 K}{\pi \delta^4}
\left( \frac{ \left|\mathbf{Y}\right| - \left|\mathbf{F}\xi\right| }{ \left|\xi\right| } \right)
\frac{\mathbf{Y}}{\left|\mathbf{Y}\right|}
\end{equation}

$A$ is a new penalty hyperparameter. It would be nice if it were
insensitive to the grid size.

---

## Influence Functions

$w(\mathbf{r})$ is another arbitrary choice!

\begin{equation}
w(x'-x)=0 \quad\text{for}\quad x' \notin H(x)
\end{equation}

I tried all of these (and more):

| Name   | $w(r)$ in ${H}(x)$          |
|-------|-------------------------------------|
| Const  | 1                                   |
| Inv    | $\delta/r$                          |
| Linear | $1-\frac{r}{\delta}$                |
| Quadr  | $\left(1-\frac{r}{\delta}\right)^2$ |
| Cubic  | $\left(1-\frac{r}{\delta}\right)^3$ |

Different meshless methods have requirements on their radial
kernel. Does PD have any?

---

## PeriFlakes

[https://github.com/afqueiruga/PeriFlakes](https://github.com/afqueiruga/PeriFlakes)

- A cleaned up implementation of the code for *Queiruga, 2017*
- Released under GPL v3
- Written using my Cornflakes/Popcorn domain specific language
- Autogenerates code for new methods, including tangent matrices
- Jupyter notebook for all the plots in this presentation

- Databases of solutions archived at:
> Alejandro Francisco Queiruga. (2018). 2D Peridynamic Displacement Fields [Data set]. Zenodo. http://doi.org/10.5281/zenodo.1284634

## Code Generation

The file `peri_kernels_pop.py` contains object-oriented implementations like this:
```python
class Oterkus(StateBased):
    def m_expr(self):
        return self.w0I * alpha_I * rxabs
    def a_expr(self):
        return self.w0I * alpha_I * rxabs**2
    def b_expr(self):
        return self.w0I * alpha_I * (xn[0] * xn[1] )**2 * rxabs**2
    def theta_expr(self):
        return 2/m[0] * self.w0I * alpha_I * (ryabs-rxabs) *  (yn.T*xn)[0,0]
    def force(self):
        a = (c_K - 2*c_G*a_integ[0]/b_denom[0])/2
        b = c_G /( 2* b_denom[0])
        A = 2 * self.w0I* alpha_I  * ( a*2/m[0]*(yn.T*xn)[0,0]*theta[0] + b*(ryabs-rxabs) )
        return A * yn
```

- Vector symbolic notation using Sympy.
- Symbolic differentiation for tangent
- Generates a C module for matrix/vector contribution.

## Code Generation

We want to exhaustively search for methods:
```python
formulations = {
    'Silling' : StateBased,
    'Oterkus2': Oterkus,
    'Fbased':Fbased,
	'Fstab_Littlewood':Fstab_Littlewood2010,
    'Fstab_Littlewood2011':Fstab_Littlewood2011,
    'Fstab_Silling':Fstab_Silling,
}
delta = i_delta[0]
weight_funcs = {
    'const' : lambda r : Piecewise( (1.0    , Le(r,delta)), (0.0,True)),
    'inv'   : lambda r : Piecewise( (delta/r, Le(r,delta)), (0.0,True)),
    'linear': lambda r : Piecewise( ((1.0-r / delta)**1.0, Le(r,delta)), (0.0,True)),
    'quadr' : lambda r : Piecewise( ((1.0-r / delta)**2.0, Le(r,delta)), (0.0,True)),
    'cubic' : lambda r : Piecewise( ((1.0-r / delta)**3.0, Le(r,delta)), (0.0,True)),
    'quarticA':lambda r: Piecewise( (1.0 - 6.0*(r/delta)**2 + 8.0*(r/delta)**3 - 3.0*(r/delta)**4, Le(r,delta)), (0.0,True))
}
for methodname,c in formulations.iteritems():
    for weightname,f in weight_funcs.iteritems():
       c(methodname+"_"+weightname, f).kernel()
```
This generates 36 different methods with only 240 lines of code!

---

## Constant Strain Modes

*Queiruga, 2017*

---

## Silling's Method Search

*Queiruga, 2017*

---

## Oterkus's Method Search

*Queiruga, 2017*

---

## F-based Method Search

*Queiruga, 2017*

---

## Convergence Orders

How quickly do those methods converge to the true solution:

<table><tr><td></td><td>isotropic</td><td>shear</td><td>uniaxial</td></tr><tr><td>Fbased</td><td>0.998518471486</td><td>0.88580762397</td><td>0.961508949756</td></tr><tr><td>Oterkus2</td><td>-0.0934425713562</td><td>0.992619961121</td><td>-0.133199906094</td></tr><tr><td>Silling</td><td>0.940215702548</td><td>-0.155718080685</td><td>-0.120152361872</td></tr></table>

- Best results out of the exhaustive search
- F-based has $O(h^1)$ convergence, which is okay but slow!
- The dilation-based models do not converge to the correct result for
  all three cases.
  - Silling's method defines equivalance by a trial dilation
  - Oterkus's method defines equivalence by a trial shear
- Cubic cone is the best influence function

---

## Surface Effects

The $x$ component of Silling's with different $m$s:

*PeriFlakes/surface_corrections.ipynb*

---

## Surface Effects

The $x$ component of Silling's with different $m$s:

*PeriFlakes/surface_corrections.ipynb*

---

## Surface Corrections

Tried different ways of handling boundary conditions:

1. No buffer zone
\begin{equation}
\mathbf{b} = \bar{\mathbf{t}}\frac{L}{A}, \quad \mathbf{u} =
\bar{\mathbf{u}},  \quad\text{on} \,\Gamma_{bound}
\end{equation}
2. A buffer of size $\delta$ that integrates the PD operator, with the
same BCs applied everywhere
\begin{equation}
\mathbf{b} = \bar{\mathbf{t}}\frac{L}{A}, \quad \mathbf{u} =
\bar{\mathbf{u}} \quad \text{in}\, \Omega_{fict}
\end{equation}
3. A buffer with a finite difference for Neumann BCs:
\begin{equation}
\frac{u^f_y-u^o_y}{\Delta y} + \nu \frac{u^+_y-u^-_y}{2 \Delta x} = \frac{t_y(1-\nu^2)}{E}
\end{equation}
4. A buffer that integrates the PD operator, and then apples the
   BC constrain.

## Surface Effects

- Boundary conditions are hard for all numerical methods, especially strong-form
- Every constant is computed numerically at each point in PeriFlakes
  - Takes into account domain edge and $\alpha$ line-of-sight
  - Described in Mandenci's 2014 book
- Le and Bobaru, 2018 tried many and found that ficticious nodes
worked

Use a finite difference stencil:
\begin{equation}
\frac{u^f_y-u^o_y}{\Delta y} + \nu \frac{u^+_y-u^-_y}{2 \Delta x} = \frac{t_y(1-\nu^2)}{E}
\end{equation}
\begin{equation}
\frac{u^f_x-u^o_x}{\Delta y} + \nu \frac{u^+_x-u^-_x}{2 \Delta x} = \frac{2 t_x (1+\nu)}{E}
\end{equation}

## The finite difference in vector notation:

Build the deformation gradient from differences,
\begin{equation}
\mathbf{F}=\frac{\left(\mathbf{u}^{f}-\mathbf{u}^{o}\right)}{\left|\mathbf{x}^{f}-\mathbf{x}^{o}\right|}\otimes\mathbf{n}+\frac{\left(\mathbf{u}^{+}-\mathbf{u}^{-}\right)}{\left|\mathbf{x}^{+}-\mathbf{x}^{-}\right|}\otimes\mathbf{e}_{t}
\end{equation}
and then apply the BC constraint by
\begin{equation}
 0=\lambda\mathrm{tr}\left[\varepsilon\right]\mathbf{n}+2G\varepsilon\mathbf{n}-\bar{\mathbf{t}}
\end{equation}
with $\varepsilon=\frac{1}{2}\left(\mathbf{F}^{T}+\mathbf{F}\right)$.

Equivalent to the stencil in Le and Bobaru, but rotates. Replace $-$
with $o$ for a corner.

---

## Surface Corrections

*PeriFlakes/surface_corrections.ipynb*

---

## F-based

- Adding the fictitious domain with same BCs still converges
 - Stencil doesn't "catch" due to instability in operator.

*PeriFlakes/surface_corrections.ipynb*

---

## Silling's dilation

- Adding the fictitious domain with same BCs still converges
- Stencil doesn't converge where it should. (Not 100% confident.)

*PeriFlakes/surface_corrections.ipynb*

---

## Oterkus's dilation

- Adding the fictitious domain with same BCs still converges
- Stencil doesn't converge where it should. (Not 100% confident.)

*PeriFlakes/surface_corrections.ipynb*

---

## Grid Convergence on Fracture

- The smoothing step is the best stabilization
- Silling's dilation is the best, but recall that it doesn't
  pass constant strain!

*Periflakes/analysis.ipynb*

---

## Conclusions

- Need to do rigorous verifications of our numerical programs.
- Evaluate programs on fundamental problems looking for expected behavior.
- Perform hyperparameter sweeps---every choice is arbitrary.

| method | Uni. | Iso. | Shear | LFM | Stable |
| - | :-: | :-: | :-: | :-: |
| Silling | ❌ | ❎ | ❌ | ✅, #1|✅|
| Oterkus | ❌ | ❌ |❎ | ❌ |✅|
| Fbased | ❎ | ❎ | ❎ | ✅ | ❌ |
| Smoothed F | ❎ | ❎ | ❎ | ✅,#2 |✅|

- Can fix the stability and improve the accuracy of the F-based
correspondence Peridynamics model.
- Silling's dilation model is the best on the crack, but it doesn't
pass all constant-strain tests.
- BCs need work, but the error is bounded by $O(h^1)$.
- The operator itself has restrictions and needs more work.

---

# Thank you

Email: afqueiruga@lbl.gov

Checkout my codes at  
[https://github.com/afqueiruga/PeriFlakes](https://github.com/afqueiruga/PeriFlakes)  
and look at the Jupyter notebook and Zenodo database.

We can churn out and test more Peridynamics models in PeriFlakes!

---

## What is our ultimate goal?

We're searching for a numerical method that:

- Passes all standard tests (continuum mechanics solutions)
- Solves the fracture problem efficiently
- Efficiently represents small changes in fractures
- Good tradeoff between accuracy and computational cost

Compare to a discretely meshed fracture solved by FEM:

- Meshed with GMSH, implemented with FEniCS
- Included in `PeriFlakes/fem/fem_crack.py`

Wall-clock time for matrix/vector assemble + solve.

---

## Error vs. Compute time

- FEM is an order of magnitude more accurate and faster.

*Periflakes/analysis.ipynb*

---

## Minimal Representable Length

- Bond breaks are discrete events at fixed locations
- Peridynamic grid spacing is the minimal length a fracture can jump
- Same as a the mesh size in node-splitting or element-deleting FEM

---

## Minimal Representable Length

FEM solves a smaller length scale faster!

*Periflakes/analysis.ipynb*