Verify Solution Satisfies a PDE: FiPy

Say I have a solution to the following implementation of the Heat Equation using FiPy:

from fipy import CellVariable, FaceVariable, TransientTerm, DiffusionTerm, Grid1D
from fipy.tools import numerix

# Set Up Grid and Time Step
nx, dx = 100, 0.01
mesh = Grid1D(dx, nx)

deltaT = 0.9 * (dx**2) / 2
steps = 1000

# Construct Solution Variable
u = CellVariable('u', mesh=mesh, rank=0)

# # Set Dirchelet Boundary Conditions
u.constrain(value=0, where=mesh.facesLeft)
u.constrain(value=0, where=mesh.facesLeft)

# # Set Initial Condition
X = mesh.cellCenters[0]
gaussInit = numerix.exp(-(X-0.5)**2 / (2*0.125**2))
u.setValue(gaussInit)

# Construct PDE
eq = TransientTerm() == DiffusionTerm(coeff=1.0)

# Solve the PDE
T = np.array([0])
Solution = u.value

time = 0
time_iter = 0
while time_iter < steps:
    time += deltaT
    eqn.solve(dt=dt, var=u)

    T = np.append(T, time)
    Solution = np.vstack((Solution, u.value))

Now that I have the solution to the PDE Solution at every time step T, I want to see if this data satisfies the PDE. How do I do that?

** Edit **

Consider the following code:

eqn = (TransientTerm(coeff=rho)
       + ConvectionTerm(coeff=v)
       == DiffusionTerm(coeff=D)
       + ImplicitSourceTerm(coeff=S1)
       + S0)

data = np.load('EquationSolution.npz')
t = data['t']
x_cell, x_face = data['x_cell'], data['x_face']
phi_cell, phi_face = data['phi_cell'], data['phi_face']
# _cell means cell-centered
# _face means face-centered
# phi_cell: solution at all times t and all locations x_cell

I want to verify that my loaded data solve the above equation eqn. I don't want to verify this by solving the PDE as that will be too costly. Rather, I want to somehow use eqn.residualVectorAndNorm to basically see if phi_cell solves eqn. How do I best go about doing this?

In general, you test against solutions you know. The FiPy documentation has lots of examples, many of which test against analytical solutions.

If you call eqn.sweep(), instead of eqn.solve(), it returns a residual, which is a measure of how well the solution satisfies the equation. I note that the documentation is wrong; it doesn't return the residual vector, it returns the L2-norm of the residual vector: , where L is the linearized discretization matrix of the PDE, x is the solution vector, and b is the right-hand-side of the equation (due to boundary conditions, sources, and "old" values). Note that the residual is evaluated before the equation is solved; after the equation is solved, the residual will be whatever the tolerance was of the linear solver, but it won't necessarily be a good solution to the nonlinear PDE. If you have non-linear coefficients, the residual should get smaller as you repeatedly sweep the equation. For your example with constant coefficients, you can sweep a second time (or call eqn.residualVectorAndNorm()).

If you don't trust our discretizations, then you shouldn't believe the residual, since they both use the same code. You can use the Method of Manufactured Solutions to see if FiPy is doing what you expect. Briefly, that involves making up a solution, analytically applying your equation to that solution, and then adding that analytical result to your equation as a source. You should get back your made-up (manufactured) solution (within discretization error; this is a way to measure discretization error).

Edit in response to Edited Edit

Starting from

I want to verify that my loaded data solve the above equation eqn. I don't want to verify this by solving the PDE as that will be too costly. Rather, I want to somehow use eqn.residualVectorAndNorm to basically see if phi_cell solves eqn. How do I best go about doing this?

You'd load old and current values you want to test by something like

phi_cell.value = data['phi_cell_old']
phi_cell.updateOld()
phi_cell.value = data['phi_cell']
res, resnorm = eqn.residualVectorAndNorm(var=phi, dt=dt)

The details will depend on exactly what you have in EquationSolution.npz.

This will discretize the terms to the semi-implicit linear matrix, and it will then explicitly solve , given your proposed value of phi_cell. It will then return both and .