Two different results from coupled eqns

[Sandip433]

I want to simulate the coupled equations of photodesorption of Oxygen molecules from palladium surface due to laser irradiation in femtosecond range . These are the equations .
Desorption.pdf
In fipy I substitute the adsorbate temperature T_ads in other equations to solve directly and I got four equations , 1 - electrons temperature, 2 - phonons temperature, 3 - Adsorbate temperature , 4 - Desorption yield.
Here is the first code

from fipy import *

Lx=1.0e-7#m
dx=1.0e-9 #m
nx=int(Lx/dx) #number of cells in mesh
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
E_a=2.15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[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)+((2/(numerix.exp((time- 0e-15)/t_0)+numerix.exp(-(time-0e-15)/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)
T_ads=CellVariable(name="Adsorbate Temperature", mesh=m, value=350.0, hasOld=True)
Theta=CellVariable(name="Change in Surface Coverage", mesh=m, value=0.5, 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)
eta_e=1e13
eta_p=1e10
k=8.6173303e-5 #eV/K
h=4.1357e-15 #eV*s
neu_ads=1e13
neu_rc=1.3e13
T_1=1/((neu_ads*(h/k))*numerix.exp((neu_ads*(h/k))/T_ads))
T_2=(T_ads*(numerix.exp((neu_ads*(h/k))/T_ads)-1))**2
T_3=eta_e*((neu_rc/neu_ads)/(numerix.exp((neu_rc*(h/k))/T_e)-1)-1/(numerix.exp((neu_ads*(h/k))/T_ads)-1))
T_4=eta_p*((neu_rc/neu_ads)/(numerix.exp((neu_rc*(h/k))/T_p)-1)-1/(numerix.exp((neu_ads*(h/k))/T_ads)-1))
neu_pfac= 1e13
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)+S
eq_p=TransientTerm(coeff=C_p, var=T_p)==ImplicitSourceTerm(coeff=G, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_p)
eq_a=TransientTerm(coeff=1, var=T_ads)== T_1*T_2*(T_3+T_4)
eq_t=TransientTerm(coeff=1, var=Theta)== Theta *neu_pfac*numerix.exp(-E_a/(k*T_ads))
eq=eq_e & eq_p & eq_a & eq_t
results_T_e = []
results_T_p = []
results_T_ads = []
results_Theta = []
t_list=[]
temperature_e_data = []
temperature_p_data = []
temperature_a_data = []
temperature_t_data = []
while time()<100e-12:
    t_list.append(float(time()))
    results_T_e.append(float(T_e.max()))
    results_T_p.append(float(T_p.max()))
    results_T_ads.append(float(T_ads.max()))
    results_Theta.append(float(Theta.max()))
    T_e.updateOld()
    T_p.updateOld()
    T_ads.updateOld()
    Theta.updateOld()
    for i in range(3):
        res=eq.sweep(dt=0.5e-15)
    time.setValue(time() + 0.5e-15)

in this case I see that there is no change in desorption yield. So I put these h= 0.5 - Theta , in the equations and I can see the change.

from fipy import *

Lx=1.0e-7#m
dx=1.0e-9 #m
nx=int(Lx/dx) #number of cells in mesh
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
E_a=2.15
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[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)+((2/(numerix.exp((time- 0e-15)/t_0)+numerix.exp(-(time-0e-15)/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)
T_ads=CellVariable(name="Adsorbate Temperature", mesh=m, value=350.0, hasOld=True)
Theta=CellVariable(name="Change in Surface Coverage", mesh=m, value=0.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)
eta_e=1e13
eta_p=1e10
k=8.6173303e-5 #eV/K
h=4.1357e-15 #eV*s
neu_ads=1e13
neu_rc=1.3e13
T_1=1/((neu_ads*(h/k))*numerix.exp((neu_ads*(h/k))/T_ads))
T_2=(T_ads*(numerix.exp((neu_ads*(h/k))/T_ads)-1))**2
T_3=eta_e*((neu_rc/neu_ads)/(numerix.exp((neu_rc*(h/k))/T_e)-1)-1/(numerix.exp((neu_ads*(h/k))/T_ads)-1))
T_4=eta_p*((neu_rc/neu_ads)/(numerix.exp((neu_rc*(h/k))/T_p)-1)-1/(numerix.exp((neu_ads*(h/k))/T_ads)-1))
neu_pfac= 1e13
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)+S
eq_p=TransientTerm(coeff=C_p, var=T_p)==ImplicitSourceTerm(coeff=G, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_p)
eq_a=TransientTerm(coeff=1, var=T_ads)== T_1*T_2*(T_3+T_4)
eq_t=TransientTerm(coeff=1, var=Theta)== (0.5-Theta) *neu_pfac*numerix.exp(-E_a/(k*T_ads))
eq=eq_e & eq_p & eq_a & eq_t
results_T_e = []
results_T_p = []
results_T_ads = []
results_Theta = []
t_list=[]
temperature_e_data = []
temperature_p_data = []
temperature_a_data = []
temperature_t_data = []
while time()<100e-12:
    t_list.append(float(time()))
    results_T_e.append(float(T_e.max()))
    results_T_p.append(float(T_p.max()))
    results_T_ads.append(float(T_ads.max()))
    results_Theta.append(float(Theta.max()))
    T_e.updateOld()
    T_p.updateOld()
    T_ads.updateOld()
    Theta.updateOld()
    for i in range(3):
        res=eq.sweep(dt=0.5e-15)
    time.setValue(time() + 0.5e-15)

then I thought that due to very small change in desorption ,this happens. Problem arises when I want to do parametric variation with respect of neu_pfac. After calculation I see that yield is linearly dependent on neu_pfac but from the equations I can see that it can not be. I could not correct the code because I did not understand where are the errors.

[guyer]

What's the question? I don't know your equations or how they're supposed to behave.

[guyer]

I don't know if it will change anything, but you'll have mildly better coupling if you convert the right-hand side of eq_a to an ImplicitSourceTerm in T_ads and the RHS of eq_t to an ImplicitSourceTerm in Theta.