CellVariable(...hasOld=0) influences results from solve()

[AnnikaSchultheiss]

I was running the 1D diffusion example. I realised, that using CellVariable(...hasOld=0) and CellVariable(...hasOld=1) give different results. From my understanding hasOld should not influence the result of solve.

from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer
from fipy.tools import numerix
import matplotlib.pyplot as plt

nx = 50
dx = 1.
steps = 100
mesh = Grid1D(nx=nx, dx=dx)
phi = CellVariable(name="solution variable",mesh=mesh,value=0., hasOld=0) #<<<<<<<<<<<<<<<<<<
D = 1.
valueLeft = 1
valueRight = 0
phi.constrain(valueRight, mesh.facesRight)
phi.constrain(valueLeft, mesh.facesLeft)
timeStepDuration = 0.9 * dx**2 / (2 * D)
eqX = TransientTerm() == DiffusionTerm(coeff=D)
for step in range(steps):
    eqX.solve(var=phi,dt=timeStepDuration)
plt.plot(mesh.x.value, phi.value)
plt.show()

CellVariable(...hasOld=0), which is the correct solution according to the analytical solution given in 1D diffusion example:

From my understanding hasOld should not influence the result of solve. But CellVariable(...hasOld=1):

Note: In the documentation of CellVariable hasOld is not described. Version: 3.4.4

[wd15]

I think it will influence the result. If hasOld=1 then the "old" value is never updated so the loop just continues to start at the old value. Hence why the solution doesn't evolve. If you set hasOld=1 and then include phi.updateOld() in the loop then things should work.

[raviapatel]

thanks there was a bit of a confusion on the updateOld. We thought it was more to save the result of previous timestep but for solving in implicit form ct is the current input variable I am also very curious in how fipy handles expressions. This is not clearly documented for instance if the coefficient is a function of primary variable or an instance where you add additional sink/source term like (phis-phis.old)/dt in equation above where phis is another indepenedent variable how does FiPy handle this if phis is changing in terms of assmbling of LHS and RHS matrix. It would be great to have some clarity on this

[guyer]

TransientTerm(coeff=rho, var=phi) represents $(\rho\phi - \rho^\text{old}\phi^\text{old}) V_P / \Delta t$. If you do not specify hasOld=True, then FiPy uses the current value of phi when it puts $\rho^\text{old}\phi^\text{old} V_P / \Delta t$ on the RHS vector and the solution evolves with each step. When you do specify hasOld=True, then phi.old is used and phi.old only changes when you explicitly call phi.updateOld(); if you don't call phi.updateOld(), then it just does the same solution over and over.

Have you seen the documentation of the general form of an equation and the descretization of the terms in that equation? Any term involving what you refer to as the primary term, e.g. phi in BlahBlahTerm(var=phi) or eq.solve(var=phi), will be rendered as described in the discretization section, with the coefficients of anything implicit in phi on the LHS matrix and everything else on the RHS vector, including phi.old. In the case of an expression like (phis-phis.old)/dt, then both phis and phis.old will be explicit and put on the RHS vector. If you use coupled equations, however, you can add TransientTerm(var=phis) to your equation and the coefficient of phis will go to the appropriate block of the LHS matrix and phis.old will go to the appropriate subvector of the RHS. More generally, in an equation for phi, adding + k * phis will be explicit (and go to the RHS) and adding + ImplicitSourceTerm(coeff=k, var=phis) will be implicit and go to the LHS matrix (if it can; it depends on the sign relative to what's already on the matrix, but FiPy handles it automatically).

[raviapatel]

Great thanks for the explaination this makes things more clear. Just a short follow up if " + k * phis will be explicit (and go to the RHS) " then if phis updates RHS matrix would automatically get updated ? Second is if the coefficient are not values but function calls does fipy call this function on every time iteration?

[guyer]

The current value of phis will be used to build the RHS vector each time you call solve() or sweep().