Connected to this question here on Computational Science, I've posted a follow-up question on how to solve a PDE using an implicit scheme like Crank-Nicholson in general in this question on SO.
But I guess it is better to ask here it is correct to solve a PDE with an implicit step in this way in general and ask the more specific coding issues on SO. So I'll try to split up the question I asked on SO.
So I've got a non-linear advection-diffusion(-reaction) PDE which I want to solve for one single timestep using Crank-Nicolson. In this example my 1D-PDE is the following:
$$\frac{T^{n+1} - T^n}{\Delta t} \approx \frac{1}{2} \left(\left(\alpha(T^n) \frac{T^n_{i-1} - 2T^n_{i} + T^n_{i+1}}{\Delta x^2} + v(T^n) \frac{T^n_{i\pm1} - T^n_{i}}{\Delta x}\right) + \left( \alpha(T^{n+1}) \frac{T^{n+1}_{i-1} - 2T^{n+1}_{i} + T^{n+1}_{i+1}}{\Delta x^2} + v(T^{n+1}) \frac{T^{n+1}_{i\pm1} - T^{n+1}_{i}}{\Delta x}\right)\right)$$
Where $\alpha$ is the temperature dependent diffusion coefficient and $v$ is the temperature dependent advection speed. The advection is discretized using an upwind scheme, since when there is advection $|Pe|\geq2$ can be assumed.
My current program code returns the derivatives of the PDE for given input temperatures, summing up the specific advection and diffusion derivatives. Thus my solver can't "see" the diffusion coefficient and the advection speed. Using a matrix notation to solve the PDE is thus not applicable.
Furthermore the derivatives and value (temperature) arrays have the same shape as the calculation grid, meaning that for the simplest case of the 1D PDE with, for example 10 nodes, the derivatives and value arrays are vectors with 10 cells (of course neglecting the additional ghost cells I use for fast stencil solving and boundary conditions).
The derivative terms in my PDE are thus in my program code replaced with the following vector derivative where each element in $derivative$ corresponds to the derivative of its respective temperature element in the value array: $$derivative(T^n) = \alpha(T^n) \frac{T^n_{i-1} - 2T^n_{i} + T^n_{i+1}}{\Delta x^2} + v(T^n) \frac{T^n_{i\pm1} - T^n_{i}}{\Delta x} $$
So in pseudo-code equation my PDE looks like this: $$\frac{T^{n+1} - T^n}{\Delta t} \approx \frac{1}{2} \left(derivative(T^{n}) + derivative(T^{n+1})\right)$$
To solve this, I reorganize this equation so that $$T^{n+1} - T^n - \frac{\Delta t}{2} \left(derivative\left(T^{n}\right) + derivative(T^{n+1})\right) = 0$$
which in python looks like
def crank_nicolson(y, yprev, h):
return (y - yprev - h / 2 * (diff(y) + diff(yprev)))
with the stepsize $h = \Delta t$, $y=T^{n+1}$, $yprev=T^{n}$ and $diff = derivative(T)$.
This can easily be fed to any root finding algorithm like the scipy.optimize.root module or using any other self-written iterative solver for example like:
from scipy import optimize
def diff(T): # simply example differential function
return np.sqrt(T)
def iterate_cranknicolson(y0, h, rtol, diff): # simple iterative solver
y = y0
diff_old = diff(y) # save forward diff
f = np.ones_like(y)
i = 0
while np.any(f > rtol):
y_old = y
y = y0 + h / 2 * (diff_old + diff(y)) # at the first step: euler forward predictor
f = (y / y_old - 1) # simple stop criterion
i += 1
return y
# some test runs
T0 = np.arange(10, 20).astype(float) # set starting temperature
cn_hybr = optimize.root(crank_nicolson, T0, args=(T0, 1.), method='hybr', tol=1e-6)
cn_it = iterate_cranknicolson(T0, 1., 1e-6, diff)
TL,DR. The question(s) summmarized:
- Is this approach to iteratively solve an implicit discretization of a PDE correct? With "this approach" I mean especially reorganizing the PDE to $T^{n+1} - T^n - \frac{\Delta t}{2} \left(derivative\left(T^{n}\right) + derivative(T^{n+1})\right) = 0$ and finding the root of this equation.
- Is there any better way to do it without being able to use matrix notation?
- Is the derivative of $T^n$ the first column of the Jacobian matrix, so that I could use a modified Newtonian iteration, as proposed here and does this work without having a matrix notation at all?
- (not really important, thus secondary: Is it ok to use upwind scheme for advection in Crank-Nicolson or should I stick to central differences?)
Thanks alot in advance! Please tell me if I should refine parts of my question and/or add more detail.
