Dear Daniel, Jonathan,
Thank you both very much for your helpful replies. Krishna and I can now access
the 'L' discretisation matrix in Lx = b, as desired, using the code example at
the bottom of this e-mail.
We have some follow-up questions, though:
1) Applying numerix.array() to 'L', when 'L' is of type
'fipy.matrices.scipyMatrix._ScipyMeshMatrix', creates a zero-dimensional
ndarray, with no shape. This isn't what we expected because L has diagonal
numerical values & '---' where its sparse "entries" are.
Our goal is to obtain 'L' using your suggested method and then to convert it
into the SciPy sparse.csc_matrix format for further processing. The input to
SciPy's csc_matrix function must be a rank-2 ndarray, but (reasonably enough!)
this fails when we pass csc_matrix a zero-dimensional 'L' matrix.
2) We see from the 2009 paper that it's a three-point stencil used for the
generation of the discretisation matrix in a first order scheme. What stencil
is used for 2nd order schemes?
3) How do we implement a higher (e.g. 8th & 12th order central-difference)
order schemes in FiPy?
Thanks again.
With best regards,
- Ian & Krishna
Modified diffusionmesh20x20.py example:
from __future__ import division
from fipy import *
from scipy.sparse import csc_matrix
SparseMatrix = LinearLUSolver()._matrixClass
nx = 3
ny = nx
dx = 1.
dy = dx
L = dx * nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = CellVariable(name = "solution variable",
mesh = mesh,
value = 1.)
D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D)
valueTopLeft = -10
valueBottomRight = 20
X, Y = mesh.faceCenters
facesTopLeft = ((mesh.facesLeft & (Y > L / 2))
| (mesh.facesTop & (X < L / 2)))
facesBottomRight = ((mesh.facesRight & (Y < L / 2))
| (mesh.facesBottom & (X > L / 2)))
phi.constrain(valueTopLeft, facesTopLeft)
phi.constrain(valueBottomRight, facesBottomRight)
v, L, b = DiffusionTerm(D)._buildMatrix(phi, SparseMatrix)
print(type(L))
print(L)
print("\n")
numerixed_L = numerix.array(L)
print(type(numerixed_L))
print(numerixed_L.shape)
print(numerixed_L.ndim)
print(numerixed_L)
-----Original Message-----
From: [email protected] [mailto:[email protected]] On Behalf Of Guyer,
Jonathan E. Dr. (Fed)
Sent: 12 January 2017 16:05
To: FIPY <[email protected]>
Subject: Re: Question on accessing internal matrices of the system being solved
By "display", I mean show a color map. It's not quantitative, but it can be
useful to understand the stencil being generated.
> On Jan 12, 2017, at 10:53 AM, Guyer, Jonathan E. Dr. (Fed)
> <[email protected]> wrote:
>
> In addition to this, if you define an environment variable
> FIPY_DISPLAY_MATRIX, then FiPy will display the matrix and RHS vector. If
> FIPY_DISPLAY_MATRIX contains "terms", then the contribution of each term will
> be displayed before the aggregate matrix and vector are displayed.
> If FIPY_DISPLAY_MATRIX contains "print", then the numerical values of the
> matrix and vector will be printed to the shell. This is overwhelmingly
> useless for any but very tiny meshes.
>
>> On Jan 12, 2017, at 10:27 AM, Daniel Wheeler <[email protected]>
>> wrote:
>>
>> Hi Krishna,
>>
>> Yes you can do this, see
>>
>> https://github.com/usnistgov/fipy/blob/develop/fipy/terms/term.py#L33
>> 2
>>
>> The following should work,
>>
>> equation.cacheMatrix()
>> equation.solve(...)
>> the_matrix = equation.matrix
>>
>> The cacheMatrix is require because otherwise the reference to the
>> matrix is lost.
>>
>> I'm not exactly certain on the matrix format, but you should be able
>> to do
>>
>> np.arrray(the_matrix)
>>
>> to make it into a 2D numpy array.
>>
>> Hope that helps.
>>
>> On Wed, Jan 11, 2017 at 3:39 PM, Gopalakrishnan, Krishnakumar
>> <[email protected]> wrote:
>>> Hi,
>>>
>>> I was wondering if there is a way to access the internal matrices of
>>> the system of equations formulated by Fipy , given the coefficient
>>> terms of the general conservation equation.
>>>
>>> Let's say I have a simple diffusion equation (eg. Heat equation). I
>>> can hand-code the differentiation matrix and BC vector for a
>>> constant dx using the 2nd order central difference formula.
>>>
>>> Is there a way to see the matrix/vector system of the discretised
>>> system that FiPy has internally formed before we call solve() or
>>> sweep() ? I am still learning basics of discretisarion schemes and
>>> it will be helpful to visualise how FiPy internally forms the
>>> matrices, and how they compare to the hand-derived matrices/vectors
>>>
>>> Regards,
>>>
>>> Krishna
>>>
>>>
>>> _______________________________________________
>>> fipy mailing list
>>> [email protected]
>>> http://www.ctcms.nist.gov/fipy
>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>>
>>
>>
>>
>> --
>> Daniel Wheeler
>> _______________________________________________
>> fipy mailing list
>> [email protected]
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
