# examples.phase.impingement.mesh40x1ΒΆ

Solve for the impingement of two grains in one dimension.

In this example we solve a coupled phase and orientation equation on a one dimensional grid. This is another aspect of the model of Warren, Kobayashi, Lobkovsky and Carter [13]

```
>>> from fipy import CellVariable, ModularVariable, Grid1D, TransientTerm, DiffusionTerm, ExplicitDiffusionTerm, ImplicitSourceTerm, GeneralSolver, Viewer
>>> from fipy.tools import numerix
```

```
>>> nx = 40
>>> Lx = 2.5 * nx / 100.
>>> dx = Lx / nx
>>> mesh = Grid1D(dx=dx, nx=nx)
```

This problem simulates the wet boundary that forms between grains of different orientations. The phase equation is given by

where

and the orientation equation is given by

where

The initial conditions for this problem are set such that for and

Here the phase and orientation equations are solved with an explicit and implicit technique respectively.

The parameters for these equations are

```
>>> timeStepDuration = 0.02
>>> phaseTransientCoeff = 0.1
>>> thetaSmallValue = 1e-6
>>> beta = 1e5
>>> mu = 1e3
>>> thetaTransientCoeff = 0.01
>>> gamma= 1e3
>>> epsilon = 0.008
>>> s = 0.01
>>> alpha = 0.015
```

The system is held isothermal at

```
>>> temperature = 1.
```

and is initially solid everywhere

```
>>> phase = CellVariable(
... name='phase field',
... mesh=mesh,
... value=1.
... )
```

Because `theta`

is an -valued variable (i.e. it maps to the circle) and thus
intrinsically has -periodicity,
we must use `ModularVariable`

instead of a `CellVariable`

. A
`ModularVariable`

confines `theta`

to
by adding or subtracting where
necessary and by defining a new
subtraction operator between two angles.

```
>>> theta = ModularVariable(
... name='theta',
... mesh=mesh,
... value=1.,
... hasOld=1
... )
```

The left and right halves of the domain are given different orientations.

```
>>> theta.setValue(0., where=mesh.cellCenters[0] > Lx / 2.)
```

The `phase`

equation is built in the following way.

```
>>> mPhiVar = phase - 0.5 + temperature * phase * (1 - phase)
```

The source term is linearized in the manner demonstrated in
`examples.phase.simple`

(Kobayashi, semi-implicit).

```
>>> thetaMag = theta.old.grad.mag
>>> implicitSource = mPhiVar * (phase - (mPhiVar < 0))
>>> implicitSource += (2 * s + epsilon**2 * thetaMag) * thetaMag
```

The `phase`

equation is constructed.

```
>>> phaseEq = TransientTerm(phaseTransientCoeff) \
... == ExplicitDiffusionTerm(alpha**2) \
... - ImplicitSourceTerm(implicitSource) \
... + (mPhiVar > 0) * mPhiVar * phase
```

The `theta`

equation is built in the following way. The details for
this equation are fairly involved, see J. A. Warren *et al.*. The main
detail is that a source must be added to correct for the
discretization of `theta`

on the circle.

```
>>> phaseMod = phase + ( phase < thetaSmallValue ) * thetaSmallValue
>>> phaseModSq = phaseMod * phaseMod
>>> expo = epsilon * beta * theta.grad.mag
>>> expo = (expo < 100.) * (expo - 100.) + 100.
>>> pFunc = 1. + numerix.exp(-expo) * (mu / epsilon - 1.)
```

```
>>> phaseFace = phase.arithmeticFaceValue
>>> phaseSq = phaseFace * phaseFace
>>> gradMag = theta.faceGrad.mag
>>> eps = 1. / gamma / 10.
>>> gradMag += (gradMag < eps) * eps
>>> IGamma = (gradMag > 1. / gamma) * (1 / gradMag - gamma) + gamma
>>> diffusionCoeff = phaseSq * (s * IGamma + epsilon**2)
```

The source term requires the evaluation of the face gradient without
the modular operator. `theta`

`faceGradNoMod`

evaluates the gradient without modular arithmetic.

```
>>> thetaGradDiff = theta.faceGrad - theta.faceGradNoMod
>>> sourceCoeff = (diffusionCoeff * thetaGradDiff).divergence
```

Finally the `theta`

equation can be constructed.

```
>>> thetaEq = TransientTerm(thetaTransientCoeff * phaseModSq * pFunc) == \
... DiffusionTerm(diffusionCoeff) \
... + sourceCoeff
```

If the example is run interactively, we create viewers for the phase and orientation variables.

```
>>> if __name__ == '__main__':
... phaseViewer = Viewer(vars=phase, datamin=0., datamax=1.)
... thetaProductViewer = Viewer(vars=theta,
... datamin=-numerix.pi, datamax=numerix.pi)
... phaseViewer.plot()
... thetaProductViewer.plot()
```

we iterate the solution in time, plotting as we go if running interactively,

```
>>> steps = 10
>>> from builtins import range
>>> for i in range(steps):
... theta.updateOld()
... thetaEq.solve(theta, dt = timeStepDuration)
... phaseEq.solve(phase, dt = timeStepDuration)
... if __name__ == '__main__':
... phaseViewer.plot()
... thetaProductViewer.plot()
```

The solution is compared with test data. The test data was created
with `steps = 10`

with a FORTRAN code written by Ryo Kobayashi for
phase field modeling. The following code opens the file `mesh40x1.gz`

extracts the data and compares it with the `theta`

variable.

```
>>> import os
>>> from future.utils import text_to_native_str
>>> testData = numerix.loadtxt(os.path.splitext(__file__)[0] + text_to_native_str('.gz'))
>>> testData = CellVariable(mesh=mesh, value=testData)
>>> print(theta.allclose(testData))
1
```