On 19 February 2014 02:25, Daniel Wheeler <[email protected]> wrote:
> The difficulty is of course that the
> coefficients are functions of the variables so this would have to be
> automated somehow, probably with sympy.
>
Many thanks for taking the time to respond. I've thought about it a
little more. Say your equation (steady-state) is
eqn = 0 == DiffusionTerm(phi)
Consider only a 1D mesh. It turns out that for this system the
(Jacobian * dx^2) is
|-3 1 0 0 0 ...
|1 -2 1 0 0 ...
|0 1 -2 1 0 ...
|0 0 1 -2 1 ...
|(continues until nx rows and nx columns)
And this continues all the way down into the bottom right-hand
corner such that J[nx,nx] * dx^2 is again -3. I found this by using
perturbation (forward differencing) to calculate the Jacobian
numerically. Of course, once I saw the pattern it became obvious
that I was being a bit silly, and I then set up a finite differencing
algebraic scheme to check the diagonal (and off-diagonal) entries \
then it makes perfect sense why this should be so.
Consider for CV n:
eqn[n] = ( (phi[n+1] - phi[n])/dx - (phi[n] - phi[n-1])/dx ) / dx
eqn[n] = phi[n+1]/dx^2 - 2 * phi[n]/dx^2 + phi[n-1]/dx^2
The gradient of this w.r.t. phi[n] is
d(eqn[n]) / phi[n] = -2 / dx^2
This is the Jacobian entry for eqn[n] and phi[n], i.e. the diagonal
of the Jacobian. The off-diagonal elements ("1") follow from the
scheme above. This makes it very easy to create a Jacobian
(for the simple case) a priori.
Question: why are the J elements for CV[zero] and CV[nx]
equal to "-3"?
I'm missing something very simple probably. How is the
discretization for the diffusion term set up (numerically) at the
boundaries? I did not find this yet in the FiPy manual.
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
