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 On Thursday, August 12, 2021 at 5:57:02 PM UTC+1 jonatha...@nist.gov wrote: > 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. > > On Aug 12, 2021, at 12:20 PM, Ben Lang <benl...@gmail.com> wrote: > > 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 <daniel....@gmail.com> > 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 <benl...@gmail.com> 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 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 >> > > -- > 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, send email to fipy+unsubscr...@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+unsubscr...@list.nist.gov.
