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) ~~~~ On Tue, Aug 10, 2021 at 6:16 AM Ben Lang <[email protected]> wrote: > > 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 > > > > -- > To unsubscribe from this group, send email to [email protected] > > View this message at https://list.nist.gov/fipy > --- > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. -- Daniel Wheeler -- To unsubscribe from this group, send email to [email protected] View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to [email protected].
