examples.diffusion.anisotropy¶

Solve the diffusion equation with an anisotropic diffusion coefficient.

We wish to solve the problem

$\frac{\partial \phi}{\partial t} = \partial_j \Gamma_{ij} \partial_i \phi$

on a circular domain centered at $(0, 0)$ . We can choose an anisotropy ratio of 5 such that

$\Gamma' = \begin{pmatrix} 0.2 & 0 \\ 0 & 1 \end{pmatrix}$

A new matrix is formed by rotating $\Gamma'$ such that

$R = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix}$

and

$\Gamma = R \Gamma' R^T$

In the case of a point source at $(0, 0)$ a reference solution is given by,

$\phi \left( X, Y, t \right) = Q \frac{ \exp \left( -\frac{1}{4 t} \left( \frac{ X^2 }{ \Gamma'_{00}} + \frac{ Y^2 }{ \Gamma'_{11}} \right) \right) }{ 4 \pi t \sqrt{\Gamma'_{00} \Gamma'_{11}} }$

where $\left(X, Y \right)^T = R \left(x, y \right)^T$ and $Q$ is the initial mass.

>>> from fipy import CellVariable, Gmsh2D, Viewer, TransientTerm, DiffusionTermCorrection
>>> from fipy.tools import serialComm, numerix

Import a mesh previously created using Gmsh.

>>> import os
>>> mesh = Gmsh2D(os.path.splitext(__file__)[0] + '.msh', communicator=serialComm) 

Set the center-most cell to have a value.

>>> var = CellVariable(mesh=mesh, hasOld=1) 
>>> x, y = mesh.cellCenters 
>>> var[numerix.argmin(x**2 + y**2)] = 1. 

Choose an orientation for the anisotropy.

>>> theta = numerix.pi / 4.
>>> rotationMatrix = numerix.array(((numerix.cos(theta), numerix.sin(theta)), \
...                                 (-numerix.sin(theta), numerix.cos(theta))))
>>> gamma_prime = numerix.array(((0.2, 0.), (0., 1.)))
>>> DOT = numerix.NUMERIX.dot
>>> gamma = DOT(DOT(rotationMatrix, gamma_prime), numerix.transpose(rotationMatrix))

Make the equation, viewer and solve.

>>> eqn = TransientTerm() == DiffusionTermCorrection((gamma,))

>>> if __name__ == '__main__':
...     viewer = Viewer(var, datamin=0.0, datamax=0.001)

>>> mass = float(var.cellVolumeAverage * numerix.sum(mesh.cellVolumes)) 
>>> time = 0
>>> dt=0.00025

>>> from builtins import range
>>> for i in range(20):
...     var.updateOld() 
...     res = 1.
...
...     while res > 1e-2:
...         res = eqn.sweep(var, dt=dt) 
...
...     if __name__ == '__main__':
...         viewer.plot()
...     time += dt

Compare with the analytical solution (within 5% accuracy).

>>> X, Y = numerix.dot(mesh.cellCenters, CellVariable(mesh=mesh, rank=2, value=rotationMatrix)) 
>>> solution = mass * numerix.exp(-(X**2 / gamma_prime[0][0] + Y**2 / gamma_prime[1][1]) / (4 * time)) / (4 * numerix.pi * time * numerix.sqrt(gamma_prime[0][0] * gamma_prime[1][1])) 
>>> print(max(abs((var - solution) / max(solution))) < 0.08) 
True

Last updated on Jun 27, 2023. Created using Sphinx 6.2.1.