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 On Thu, Aug 12, 2021 at 2:25 AM Daniel Wheeler <[email protected]> wrote: > 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].
