Changing sweep number changes solution

[Sandip433]

I want to solve these coupled equations

This is the code I wrote

from fipy import *
import matplotlib.pyplot as plt
from tqdm import tqdm
Lx=1.0e-7#m
dx=1.0e-11 #m
nx=int(Lx/dx) #number of cells in mesh
T=2e-12
dt=4e-15
G=8.95e17  #W/m^3/K
t_0=40e-15  #s
w_0=1.24489426e-8  #beam waist in m
J=28.6  #Fluence J/m^2
gamma1=78849  #J/m^3/K^2
gamma2=249.14  #J/m^3/K^3
K_0=72.0  #W/m/K
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
X = m.faceCenters[0]
time = Variable(0.)
S=(J/(2*t_0*w_0))*numerix.exp(-x/w_0)*((2/(numerix.exp(time/t_0)+numerix.exp(-time/t_0)))**2)
T_e=CellVariable(name="Electron Temperature", mesh=m, value=350.0, hasOld=True)
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=350.0, hasOld=True)
C_p=(2801128 + (-357785-2801128)/(1+(T_p/62.98191)**2.06271) + 278.44328*T_p)  #J/m^3/K
K_e=K_0*(T_e/T_p)
eq_e=TransientTerm(coeff=gamma1+gamma2*T_e, var=T_e)==DiffusionTerm(coeff=K_e, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_e)+ImplicitSourceTerm(coeff=G, var=T_p)+ImplicitSourceTerm(coeff=S/T_e,var=T_e)
eq_p=TransientTerm(coeff=C_p, var=T_p)==ImplicitSourceTerm(coeff=G, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_p)
eq=eq_e & eq_p 
results_T_e = []
results_T_p = []
t_list=[]
progress_bar = tqdm(total=int(T / dt), desc='Solving PDE')
while time()<T:
    t_list.append(float(time()))
    results_T_e.append(float(T_e[0]))
    results_T_p.append(float(T_p[0]))
    T_e.updateOld()
    T_p.updateOld()
    for sweep in range(1):
        res=eq.sweep(dt=dt)
    time.setValue(time() + dt)
    progress_bar.update() 
progress_bar.close() 
plt.plot(t_list,results_T_e)
plt.plot(t_list,results_T_p)

Now when the sweep number is 1 I get expected result (similar to the result reported in the literature) but when I change sweep number to 2 I get drastically different result. It seems that the effect of sweep number is divided into two sequence, odd number sequence and even number sequence. In the sweep number 1 among the odd number sweep is giving highest value and 2 among even is giving lowest value. Odd number sweeps give decreasing result and even number sweeps give increasing results. After 20-21 sweep the results converged to a stable result. Why is that happening?

[guyer]

The coefficients of your equations are very nonlinear functions of the solution variables. .sweep() or .solve() is obtaining an implicit solution to the linear discretization of your PDEs. Sweeping is used to obtain consistent solutions to the nonlinear equations.

Ideally, you should sweep until the residual no longer decreases. In your case, I find 50 to 100 sweeps are necessary for that to happen, which is a lot. You can try taking smaller time steps or you can evaluate the solutions after fewer sweeps (maybe 5 or 10) and see if the solutions are acceptable. The interplay between time steps and sweep counts are more art than science.

[Sandip433]

With only 1 sweep I get the expected and acceptable result which was mostly calculated using finite difference scheme in many literature article and I am getting same results changing timestep from 4 fs to 1 fs to 0.5 fs. Will you please guide me how to modify my code to check how the residual is decreasing with every sweep ?

[guyer]

Your code says res=eq.sweep(dt=dt). Print res.

[Sandip433]

Now I have changed the equations from

eq_e=TransientTerm(coeff=gamma1+gamma2*T_e, var=T_e)==DiffusionTerm(coeff=K_e, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_e)+ImplicitSourceTerm(coeff=G, var=T_p)+ImplicitSourceTerm(coeff=S/T_e,var=T_e)

eq_p=TransientTerm(coeff=C_p, var=T_p)==ImplicitSourceTerm(coeff=G, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_p)

to

eq_e=TransientTerm(coeff=1, var=T_e)==DiffusionTerm(coeff=K_e/C_e, var=T_e)-ImplicitSourceTerm(coeff=G/C_e, var=T_e)+ImplicitSourceTerm(coeff=G/C_e, var=T_p)+ImplicitSourceTerm(coeff=S/(C_e*T_e),var=T_e)

eq_p=TransientTerm(coeff=1, var=T_p)==ImplicitSourceTerm(coeff=G/C_p, var=T_e)-ImplicitSourceTerm(coeff=G/C_p, var=T_p)

I find that
(1) the equations give stable result in 1 sweep and the result does not change dramatically in 10,20,21 sweeps and in all sweeps I get expected result.
(2 ) reducing timestep also does not alter the result unreasonably
(3) residual also reduced to ~0.000420 after 10 sweeps
The only concern I have that is that representation of heat equation correct to solve although it produces reasonable result?

[guyer]

For constant heat capacity, that rearrangement gets done all the time, but your heat capacity is not constant, in either space or time. Maybe it doesn't matter, if sweeping doesn't change the result.

Properly,

$$ \begin{align} \frac{\partial C_e T_e}{\partial t} &= \nabla\cdot\left(\kappa_e \nabla T_e\right) - G \left(T_e - T_p\right) + S \\ C_e\frac{\partial T_e}{\partial t} + T_e\frac{\partial C_e}{\partial t}&= \nabla\cdot\left(\kappa_e \nabla T_e\right) - G \left(T_e - T_p\right) + S \\ \frac{\partial T_e}{\partial t} + \frac{T_e}{C_e}\frac{\partial C_e}{\partial t} &= \frac{1}{C_e}\nabla\cdot\left(\kappa_e \nabla T_e\right) - \frac{G}{C_e} \left(T_e - T_p\right) + \frac{S}{C_e} \\ \frac{\partial T_e}{\partial t} + \frac{T_e}{C_e}\frac{\partial C_e}{\partial t} &= \nabla\cdot\left(\frac{\kappa_e}{C_e} \nabla T_e\right) - \kappa_e \nabla T_e \cdot\nabla\frac{1}{C_e} - \frac{G}{C_e} \left(T_e - T_p\right) + \frac{S}{C_e} \end{align} $$

so you're discarding $\frac{T_e}{C_e}\frac{\partial C_e}{\partial t}$ and $\kappa_e \nabla T_e \cdot\nabla\frac{1}{C_e}$. Similarly for the $T_p$ equation.

What you can do is try these equations you have, with the missing terms, using Newton's method. Certainly missing terms in the Jacobian get sorted out by Newton-Raphson iterations; I'm less sure about missing terms in the equations themselves.

I wouldn't do ImplicitSourceTerm(coeff=S/(C_e*T_e),var=T_e) at all. It's just S/C_e.