How to evaluate the integral

[Sandip433]

I am trying to solve these 3 coupled equations but I could not figure out how to evaluate (1)

$a_{1}, a_{2}$ are functions of $T_{l}$

[guyer]

• You have five equations. Which three are you trying to couple?
• What does $I^{n-1}$ denote? Is $a_n I^{n-1}$ an inner product of some kind? Or, for $n = 1$, is it a power such that $a_1 I^0 \equiv a_1$?
• Do you want to know how to integrate something over space? phi.cellVolumeAverage represents $\int_{x_L}^{x^R} \phi\,dx / \int_{x_L}^{x^R} \,dx$
• What have you tried so far?

[Sandip433]

(1)I am trying to couple equation (2), (3), (4) although equation (1) is also coupled with them. Actually equation (1) is an integral equation but in order to solve for the $T_{e}, T_{l}$ and N we need to solve equation (1) to get the value of I

(2)It is simple multiplication $a_{1} + a_{2} * I(t,x') + \theta * N(t, x')$ where $a_{1}$ is dependent on $T_{l}$

(3) I want to solve equation (1) with other 3 equations.

(4) I have written equation (2) to (4) in fipy syntax. But I could not represent equation (1) in fipy syntax

eq_e=TransientTerm(coeff=C_e, var=T_e)== DiffusionTerm(coeff=K_e, var=T_e) + DiffusionTerm(coeff=E_g * D_0, var=N) - ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_e) + ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_l) + ImplicitSourceTerm(coeff=theta * I, var=N) + alpha * I + beta * I**2 

eq_l=TransientTerm(coeff=C_l, var=T_l)==DiffusionTerm(coeff=K_l, var=T_l) + ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_e)-ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_l)

eq_n=TransientTerm(coeff=1, var=N)== DiffusionTerm(coeff=D_0, var=N) + ImplicitSourceTerm(coeff=(- gamma * N ** 2 + delta), var= N)+ (a_1/(h_bar * omega)) * I+(a_2/(2 * h_bar * omega)) * I**2

eq=eq_e & eq_l & eq_n

[guyer]

Thank you for your explanations.

It's probably possible to make this work with .cellVolumeAverage, but I suspect it would be very inefficient. Instead, for a function $f(t, x)$, the integral $\int_0^x f(t, x'),dx'$ can be calculated with:

from scipy import integrate

f_int = fp.CellVariable(mesh=mesh, value=integrate.cumulative_trapezoid(f, mesh.x, initial=0.))

This will need to be recalculated each time you sweep.

Here's a prototype Variable that will automatically update the integral when the integrand changes:

from scipy import integrate

class Integration1DVariable(fp.CellVariable):
    def __init__(self, var, name=''):
        r"""
        Produces a cummulative integral of the supplied function.

        Parameters
        ----------
        var : array_like or ~fipy.variables.variable.Variable
            The collection of values to integrate.
        """
        fp.CellVariable.__init__(self, mesh=var.mesh, name=name)
        self.var = self._requires(var)

    def _calcValue(self):
        return integrate.cumulative_trapezoid(self.var, self.var.mesh.x, initial=0.)

This only works in 1D and it assumes that the coordinates of the mesh are monotonic. They will be for a Grid1D, but FiPy does not generally make any assumptions about the order of cells in space.

[Sandip433]

So if I want to implement this idea, should it be done like this way ?

import fipy as fp
from fipy.solvers import LinearPCGSolver
from scipy import integrate
import numpy as np

Lx=1.0e-5 #m
dx=1.0e-8 #m
nx=int(Lx/dx) 
mesh=Grid1D(nx=nx, dx=dx)
x = mesh.cellCenters[0]
time = Variable(0.)
T=1e-12
dt=4e-15

N=CellVariable(name="Free Electron", mesh=mesh, value=1e22, hasOld=True)
T_e=CellVariable(name="Electron Temperature", mesh=mesh, value=300.0, hasOld=True)
T_l=CellVariable(name="Phonon Temperature", mesh=mesh, value=300.0, hasOld=True)

class Integration1DVariable (fp.CellVariable):
    def __init__(self, var, N, theta, a_1, a_2, name=''):
        fp.CellVariable.__init__(self, mesh=var.mesh, name=name)
        self.var = self._requires(var)
        self.N = self._requires(N)
        self.theta = theta
        self.a_1 = a_1
        self.a_2 = a_2

    def _calcValue(self):
        x = self.var.mesh.x
        integrand = self.a_1 + self.a_2 * self.var + self.theta * self.N
        integral_x = integrate.cumulative_trapezoid(integrand, x, initial=0.)
        return integral_x

integral_var = Integration1DVariable(I, N, theta, a_1, a_2, name="Integral Term")

I_0 = 1e4
I = fp.CellVariable(name="Integral", mesh=m, hasOld=True)
I.setValue = I_0 * np.exp(-integral_var.value)

eq_e=TransientTerm(coeff=C_e, var=T_e)== DiffusionTerm(coeff=K_e, var=T_e) + DiffusionTerm(coeff=E_g * D_0, var=N) - ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_e) + ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_l) + ImplicitSourceTerm(coeff=theta * I, var=N) + alpha * I + beta * I **2 
eq_l=TransientTerm(coeff=C_l, var=T_l)==DiffusionTerm(coeff=K_l, var=T_l) + ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_e)-ImplicitSourceTerm(coeff=3 * K_b * N/ tao_e, var=T_l)
eq_n=TransientTerm(coeff=1, var=N)== DiffusionTerm(coeff=D_0, var=N) + ImplicitSourceTerm(coeff=(- gamma * N **2 + delta), var= N)+ (a_1/(h_bar * omega)) * I+(a_2/(2 * h_bar * omega)) * I **2
eq=eq_e & eq_l & eq_n 

t_list=[]
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
while time()<T:
    T_e.updateOld()
    T_l.updateOld()
    N.updateOld()
    t_list.append(float(time()))
    for sweep in range(1):
        residual=eq.sweep(dt=dt,solver=solver)
        I.value = (I.old + (I + I.old) * dt / 2).value
    time.setValue(time() + dt)

[guyer]

Not quite. I would leave the definition of Integration1DVariable as I had it:

class Integration1DVariable(fp.CellVariable):
    def __init__(self, var, name=''):
        r"""
        Produces a cummulative integral of the supplied function.

        Parameters
        ----------
        var : array_like or ~fipy.variables.variable.Variable
            The collection of values to integrate.
        """
        fp.CellVariable.__init__(self, mesh=var.mesh, name=name)
        self.var = self._requires(var)

    def _calcValue(self):
        return integrate.cumulative_trapezoid(self.var, self.var.mesh.x, initial=0.)

and then write:

integral_var = Integration1DVariable(var=a_1 + a_2 * I + theta * N, name="Integral Term")
I_new = I_0 * fp.numerix.exp(-integral_var)

Then, in your sweep loop:

    for sweep in range(1):
        residual=eq.sweep(dt=dt,solver=solver)
        I.value = I_new.value

I'm not sure what the intent of this is? Are you trying to do some sort of damped update?

 I.value = (I.old + (I + I.old) * dt / 2).value

Regardless, you never call I.updateOld(), so I.old will always be its initial value, I think.

[Sandip433]

I updated the code according to your instructions and the code runs without any error but output result is not converging well . Is there any problem with sweep loop ?

from fipy import *
from numpy import *
from fipy.solvers import LinearPCGSolver
from scipy import integrate
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(['science','no-latex','notebook'])
from tqdm import tqdm
Lx=1.0e-5 #m
dx=1.0e-8 #m
nx=int(Lx/dx) 
mesh=Grid1D(nx=nx, dx=dx)
x = mesh.cellCenters[0]
time = Variable(0.)
T=1e-12
dt=4e-15
N=CellVariable(name="Free Electron", mesh=mesh, value=1e22, hasOld=True)
T_e=CellVariable(name="Electron Temperature", mesh=mesh, value=300.0, hasOld=True)
T_l=CellVariable(name="Phonon Temperature", mesh=mesh, value=300.0, hasOld=True)
I = CellVariable(name="Integral", mesh=mesh,  hasOld=True)

lamda = 1064e-9
a_1 = -5895 + 62.26 * T_l - 0.2309 * T_l ** 2 + 3.186e-4 * T_l ** 3 - 9.967e-8 * T_l ** 4
a_2=2e-11
h_bar=1.054e-34 
omega =299792458 / lamda
K_b = 1.380649e-23 
theta = 5.1e-22 * T_l / 300
gamma = 1/(1/(3.8e-43 * N ** 2) + 6e-12)
E_g = 1.86e-19
delta = 1e11 * (3 * 1.602e-19 * N * K_b * T_e - E_g) 
C_l = 2.2368e6
K_l = 1.521e5 * T_l ** (-1.226) 
tao_e = 240e-15 * (1 + (N/6e26) ** 2)
K_e = -0.556 + 7.13e-3 * T_e
D_0 = 2.98e-7 * T_e
C_e = 3 * N * K_b

class Integration1DVariable(CellVariable):
    def __init__(self, var, name=''):
        CellVariable.__init__(self, mesh=var.mesh, name=name)
        self.var = self._requires(var)
    def _calcValue(self):
        return integrate.cumulative_trapezoid(self.var, self.var.mesh.x, initial=0.)
integral_var = Integration1DVariable(var=a_1 + a_2 * I + theta * N, name="Integral Term")
I_0 = 1e4
I_new = I_0 * exp(-integral_var)

eq_e=TransientTerm(coeff=C_e, var=T_e)== DiffusionTerm(coeff=K_e, var=T_e) + DiffusionTerm(coeff=E_g * D_0, var=N) - ImplicitSourceTerm(coeff=3 * K_b * N / tao_e, var=T_e) + ImplicitSourceTerm(coeff=3 * K_b * N / tao_e, var=T_l) + ImplicitSourceTerm(coeff=theta * I, var=N) + a_1 * I + a_2 * I ** 2 
eq_l=TransientTerm(coeff=C_l, var=T_l)==DiffusionTerm(coeff=K_l, var=T_l) + ImplicitSourceTerm(coeff=3 * K_b * N / tao_e , var=T_e)-ImplicitSourceTerm(coeff=3 * K_b * N / tao_e, var=T_l)
eq_n=TransientTerm(coeff=1, var=N)== DiffusionTerm(coeff=D_0, var=N) + ImplicitSourceTerm(coeff=(- gamma  + delta), var= N)+ (a_1/(h_bar * omega)) * I+(a_2/(2 * h_bar * omega)) * I ** 2
eq=eq_e & eq_l & eq_n 
t_list=[]
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=3000)
progress_bar = tqdm(total=int(T/dt), desc='Solving PDE')
while time()<T:
    T_e.updateOld()
    T_l.updateOld()
    N.updateOld()
    t_list.append(float(time()))
    for sweep in range(10):
        residual=eq.sweep(dt=dt,solver=solver)
        I.value = I_new.value
    time.setValue(time() + dt)
    progress_bar.update()
progress_bar.close()
plt.plot(x,N.value)

[guyer]

but output result is not converging well

Can you tell me how you concluding that? From what I can tell, by looking at the residual, your initial condition is converged. It's not evolving, but that's because everything is essentially constant.

If I increase the value of a_2 to, e.g, 2e4, that gives larger variation in I, and then I start to see things happening with T_e.