Off-Diagonal elements of diffusion tensor

[Ben Lang]

Dear Fipy Users,

Summary: The matrix associated with DiffusionTerm([[[0, 1], [1, 0]]]) appears to be all-zeros. Which is not the behaviour I would expect.

Wider context: I am ultimately trying to solve a 2D problem with terms that go as the third derivative of my variable, W, with respect to position. This can be achieved using three additional variables representing the dx^2, dy^2 and dxdy derivatives of my real field variable, and then including convection terms in these derivative terms in the evolution. So one of the quantities I want is the second "mixed partial derivative" of my variable:

d^2 W / dxdy

As explained in the manual / FAQ I am able to specify the diffusion tensor. For example Including  DiffusionTerm([[[1, 0], [0, 0]]]) and DiffusionTerm([[[0, 0], [0, 1]]]) to gather the dx^2 and dy^2 derivatives works as expected.

However I am having trouble with, DiffusionTerm([[[0, 1], [1, 0]]]). I expect this to give 2 dxdy (twice the second partial derivative with respect to x and y). However it seems to do nothing at all, evidenced by the complete emptiness of the matrix associated with it.

Example code:

nx = 20

ny = 20
dx = 1E-2

dy = 1E-2

mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

W = CellVariable(name=r'$W$', mesh=mesh)

eq = DiffusionTerm( coeff= [[[0, 1], [1, 0]]], var=W)

eq.cacheMatrix()

eq.solve()

>>> Factor is exactly singular # Ignore this, we only solved it so it would get the matrix.

matrix = eq._matrix.matrix

matrix 

>>> <400x400 sparse matrix of type '<class 'numpy.float64'>'

with 0 stored elements in Compressed Sparse Row format>

The matrix is all zeros, implying the term is doing nothing.

Perhaps I am misunderstanding the role of the off-diagonals in the diffusion tensor? Or some other error (possibly a Cellvariable vs. Facevariable issue). Any help would be greatly appreciated.

Thank you very much,

All the best,

Ben

[Daniel Wheeler]

Hi Ben,

I looked into this and it seems that the anisotropic off diagonal
values are added as a source term rather than to the matrix. See
these two lines in abstractDiffusionTerm.py.

- https://github.com/usnistgov/fipy/blob/master/fipy/terms/abstractDiffusionTerm.py#L118

and

- https://github.com/usnistgov/fipy/blob/master/fipy/terms/abstractDiffusionTerm.py#L412

I can't remember fully why it is implemented this way. It is likely a
simpler implementation than using a fully implicit method.

The following code is your code reworked to show something on the RHS
vector. It has large values in the RHS vector. If the diffusion term
coeff is set to zero then there are only tiny values on the RHS vector
demonstrating that at least something is changing due to the non-zero
off diagonal diffusion coefficient values.

~~~~
from fipy import Variable, FaceVariable, CellVariable, Grid1D,
ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer,
AdvectionTerm, PowerLawConvectionTerm, VanLeerConvectionTerm, Grid2D
from fipy.tools import numerix
import numpy as np

nx = 3
ny = 3
dx = 1.
dy = 1.

mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
W = CellVariable(name=r'$W$', mesh=mesh)

W[:] = np.random.random(9)
print('W:',W)
coeff = FaceVariable(mesh=mesh, rank=2, value=0.0)
coeff[0, 0] = 0.
coeff[0, 1] = 1.
coeff[1, 1] = 0.
coeff[1, 0] = 1.

eq = (TransientTerm(1e-9) == DiffusionTerm([[[0, 1], [1, 0]]]))
eq.cacheMatrix()
eq.cacheRHSvector()
eq.solve(W, dt=1.0)

matrix = eq._matrix.matrix

print(matrix)
print(eq._RHSvector)
~~~~

> --
> To unsubscribe from this group, send email to fipy+uns...@list.nist.gov
>
> View this message at https://list.nist.gov/fipy
> ---
> To unsubscribe from this group and stop receiving emails from it, send an email to fipy+uns...@list.nist.gov.



--
Daniel Wheeler

[Ben Lang]

Dear Daniel Wheeler,

Thank you very much for the very helpful response. It being included as a source term in the RHS vector explains the confusions I have been having. 

Thank you again for you help,

Ben

[Guyer, Jonathan E. Dr. (Fed)]

If you show us your equations and where they came from, we may be able to propose an alternative form, e.g., a pair of 2nd order equations.

I guess one of those equations will be purely hyperbolic, which FiPy isn’t good at.

[Ben Lang]

Dear  Jonathan,

For my purposes I think its all good. When I solve the equations the inclusion of the off-diagonals as source terms should work (I can't say for sure as I am still testing my model). There is indeed an alternative way of framing the problem where the W_xy mixed derivative terms do not feature. That alternative is probably where I am going.

For anyone interested in the problem and how it can be re-arranged: 

The equation I am approximating is the Fokker Planck equation for the Winger function W(x,y). For  H(x, y) some (fixed) potential energy function.  The equation has third derivative terms that look like

dW/dt =  convection/diffusion stuff.....  +  H_xxx W_yyy - 3 H_xxy W_xyy + 3 H_xyy W_xxy - H_yyy W_xxx           (with the foot scripts indicating partial derivatives).

My plan to solve this using Fipy had been to define three new quantities,   W_xx,  W_yy  and W_xy.   The intent had been to use: 

eq_xx = ( ImplicitSourceTerm(coeff=1, var= W_xx ) == DiffusionTerm( coeff= [[[1, 0], [0, 0]]], var=W) ) 

eq_pp = ( ImplicitSourceTerm(coeff=1, var= W_pp ) == DiffusionTerm( coeff= [[[0, 0], [0, 1]]], var=W) )

 and      

eq_xp = ( ImplicitSourceTerm(coeff=1, var= W_xp ) == DiffusionTerm( coeff= [[[0, 0.5], [0.5, 0]]], var=W) )

Then include new convection terms in the evolution of W  that go as   Grad  .   of   D_xx W_yy  ,   - 2  D_xy W_xy   and  D_yy W_xx    ,  where D = [ - H_y, H_x ], (so D_xx = [ - H_xxy, H_xxx]).  The reason for having three new variables (and three new convection terms) was that applying the chain rule (with the Grad) on this particular set of three results in the expression above without any source terms appearing.  However if the  W_xy  terms are troublesome I can instead use only  W_xx and W_yy  terms (putting some 3's in the Drift coefficients) along with a source term to cancel the chain rule issues. This might be a cleaner way of doing it anyway as it only introduces two new variables.

Thank you very much,

All the best,

Ben

[Guyer, Jonathan E. Dr. (Fed)]

Sounds like you’re on top of things. Let us know if we can be of further help.

- Jon