Solving coupled two temperature model with FiPy

[karnapravin]

I was trying to get a code for two temperature model with FiPy. I am facing some issues such as updating the heat capacity with temperature as well as getting the final plot. I have posted the equations as well as the code which I have currently.

Thank you for your help.

Code

from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer

from NTMpy import NTMpy as ntm
from matplotlib import pyplot as plt
import numpy as np
import numericalunits as u
u.reset_units('SI')

import math
import scipy
import numpy as np
from scipy.linalg import expm, sinm, cosm

m = Grid1D(nx=100, Lx=500*1e-9)

v0 = CellVariable(mesh=m, hasOld=True, value=300)
v1 = CellVariable(mesh=m, hasOld=True, value=300)

v0.constrain(0, m.facesLeft)
v0.constrain(1, m.facesRight)

v1.constrain(1, m.facesLeft)
v1.constrain(0, m.facesRight)

#copper
length_cu   = 500*1e-9 #Length of the Material
n_cu        = 1.2878+2.2596j; #Complex refractive index
k_el_cu     = 391;#Heat conductivity
k_ph_cu   = 7;#Heat conductivity
rho_cu      = 1e3*8.96;#Density
Te=300;
Tl=300
C_el_cu     =  96.8 *Te #Electron heat capacity
C_lat_cu    = 3.62e6; #Lattice heat capacity
G_cu        = 1e17;        #Lattice electron coupling constant

x=np.linspace(0.1e-11,500e-9,100)
y=np.linspace(0.1e-11,500e-9,100)
z=np.linspace(0.1e-11,500e-9,100)
tt=np.linspace(0.1e-20,10e-12,100)
F=170
S_a=0.49
delta=15e-9
FWHM=650e-9
tp=100e-15
S=(S_a/delta)*F*np.exp((-z/delta)-((4*math.log(2))*(((x*x)+(y*y))/(FWHM*FWHM))))*(1/tp)*math.sqrt((4*math.log(2)/math.pi))*np.exp((-4*math.log(2))*(((tt-2*tp)*(tt-(2*tp))/(tp*tp))))

v0.value = 300
v1.value = 300

eqn0 = TransientTerm(var=v0) == DiffusionTerm(k_el_cu/C_el_cu, var=v0) - G_cu*(Te-Tl)+S
eqn1 = TransientTerm(var=v1) == DiffusionTerm(k_ph_cu /C_lat_cu, var=v1) + G_cu*(Te-Tl)

eqn = eqn0 & eqn1

vi = Viewer((v0, v1))

from builtins import range
for t in range(1):
    v0.updateOld()
    v1.updateOld()
    eqn.solve(dt=1.e-12)
    vi.plot()

while running i get these errors as well

[guyer]

What is your question, exactly? What errors?

[karnapravin]

Sorry I have updated the previous code, but the final plots were not coming and the plot after

vi = Viewer((v0, v1))
was a straight line.

from builtins import range
for t in range(1):
v0.updateOld()
v1.updateOld()
eqn.solve(dt=1.e-12)
vi.plot()

This section was not running probably due to dimension errors.

My question is since the equations that I posted is a third order PDE. i tried to code it by looking at FiPy example of diffusion problem. Could you please look at the code and suggest me the problem in the code. It is really important for me.
Thank you for your time

[guyer]

• These are second order PDEs.
• What is your question? Please be specific. Right now, your question sounds like "do it for me".

[karnapravin]

Sorry I did not mean like that

here is my updated code and the error.

Code

from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer
from fipy import *
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(['science','no-latex','notebook'])
import math
import scipy
import numpy as np
from scipy.linalg import expm, sinm, cosm
#copper
lx   = 500e-9 #Length of the Material
dx=1e-11
nx=lx/dx
n_cu        = 1.2878+2.2596j; #Complex refractive index
ke     = 391;#Heat conductivity
kp   = 7;#Heat conductivity
rho_cu      = 1e3*8.96;#Density
Te=300;
Tl=300
Ce     =  96.8 *Te #Electron heat capacity
Cl    = 3.62e6; #Lattice heat capacity
G        = 1e17;        #Lattice electron coupling constant

m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
time = Variable(0.)
x=np.linspace(0.1e-11,500e-9,100)
y=np.linspace(0.1e-11,500e-9,100)
z=np.linspace(0.1e-11,500e-9,100)
tt=np.linspace(0.1e-20,10e-12,100)
F=1.7
S_a=0.49
delta=15e-9
w0=650e-9
tp=100e-15
S=(S_a/delta)*F*np.exp((-z/delta)-((4*math.log(2))*(((x*x)+(y*y))/(w0*w0))))*(1/tp)*math.sqrt((4*math.log(2)/math.pi))*np.exp((-4*math.log(2))*(((tt-2*tp)*(tt-(2*tp))/(tp*tp))))
T_e=CellVariable(name="Electron Temperature", mesh=m, value=Te, hasOld=True)
T_p=CellVariable(name="Phonon Temperature", mesh=m, value=Tl, hasOld=True)
eqn0 = TransientTerm(coeff=Ce*T_e, var=T_e) == DiffusionTerm(coeff=ke, var=T_e) - ImplicitSourceTerm(coeff=G, var=T_e)+ImplicitSourceTerm(coeff=G, var=T_p)+S
eqn1 = TransientTerm(coeff=Cl*T_p, var=T_p) == DiffusionTerm(coeff=kp, var=T_p) + ImplicitSourceTerm(coeff=G, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_p)
eqn = eqn0 & eqn1
results_T_e = []
results_T_p = []
t_list=[]
while time()<4e-12:
    t_list.append(float(time()))
    results_T_e.append(float(T_e.max()))
    results_T_p.append(float(T_p.max()))
    T_e.updateOld()
    T_p.updateOld()
    for i in range(1):
        res=eqn.sweep(dt=0.5e-15)
    time.setValue(time() + 0.5e-15)

Error

ypeError                                 Traceback (most recent call last)
Cell In[260], line 8
      6 T_p.updateOld()
      7 for i in range(1):
----> 8     res=eqn.sweep(dt=0.5e-15)
      9 time.setValue(time() + 0.5e-15)

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\term.py:219, in Term.sweep(self, var, solver, boundaryConditions, dt, underRelaxation, residualFn, cacheResidual, cacheError)
    180 def sweep(self, var=None, solver=None, boundaryConditions=(), dt=None, underRelaxation=None, residualFn=None, cacheResidual=False, cacheError=False):
    181     r"""
    182     Builds and solves the `Term`'s linear system once. This method
    183     also recalculates and returns the residual as well as applying
   (...)
    217         The residual vector :math:`\vec{r}=\mathsf{L}\vec{x} - \vec{b}`
    218     """
--> 219     solver = self._prepareLinearSystem(var=var, solver=solver, boundaryConditions=boundaryConditions, dt=dt)
    220     solver._applyUnderRelaxation(underRelaxation=underRelaxation)
    221     residual = solver._calcResidual(residualFn=residualFn)

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\term.py:124, in Term._prepareLinearSystem(self, var, solver, boundaryConditions, dt)
    121         from fipy.viewers.matplotlibViewer.matplotlibSparseMatrixViewer import MatplotlibSparseMatrixViewer
    122         Term._viewer = MatplotlibSparseMatrixViewer()
--> 124 var, matrix, RHSvector = self._buildAndAddMatrices(var,
    125                                                    self._getMatrixClass(solver, var),
    126                                                    boundaryConditions=boundaryConditions,
    127                                                    dt=dt,
    128                                                    transientGeomCoeff=self._getTransientGeomCoeff(var),
    129                                                    diffusionGeomCoeff=self._getDiffusionGeomCoeff(var),
    130                                                    buildExplicitIfOther=self._buildExplcitIfOther)
    132 self._buildCache(matrix, RHSvector)
    134 solver._storeMatrix(var=var, matrix=matrix, RHSvector=RHSvector)

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\coupledBinaryTerm.py:82, in _CoupledBinaryTerm._buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     78 for varIndex, tmpVar in enumerate(var.vars):
     80     SparseMatrix.varIndex = varIndex
---> 82     tmpVar, tmpMatrix, tmpRHSvector = uncoupledTerm._buildAndAddMatrices(tmpVar,
     83                                                                          SparseMatrix,
     84                                                                          boundaryConditions=(),
     85                                                                          dt=dt,
     86                                                                          transientGeomCoeff=uncoupledTerm._getTransientGeomCoeff(tmpVar),
     87                                                                          diffusionGeomCoeff=uncoupledTerm._getDiffusionGeomCoeff(tmpVar),
     88                                                                          buildExplicitIfOther=buildExplicitIfOther)
     90     termMatrix += tmpMatrix
     91     termRHSvector += tmpRHSvector

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\binaryTerm.py:28, in _BinaryTerm._buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     24 RHSvector = 0
     26 for term in (self.term, self.other):
---> 28     tmpVar, tmpMatrix, tmpRHSvector = term._buildAndAddMatrices(var,
     29                                                                 SparseMatrix,
     30                                                                 boundaryConditions=boundaryConditions,
     31                                                                 dt=dt,
     32                                                                 transientGeomCoeff=transientGeomCoeff,
     33                                                                 diffusionGeomCoeff=diffusionGeomCoeff,
     34                                                                 buildExplicitIfOther=buildExplicitIfOther)
     36     matrix += tmpMatrix
     37     RHSvector += tmpRHSvector

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\binaryTerm.py:28, in _BinaryTerm._buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     24 RHSvector = 0
     26 for term in (self.term, self.other):
---> 28     tmpVar, tmpMatrix, tmpRHSvector = term._buildAndAddMatrices(var,
     29                                                                 SparseMatrix,
     30                                                                 boundaryConditions=boundaryConditions,
     31                                                                 dt=dt,
     32                                                                 transientGeomCoeff=transientGeomCoeff,
     33                                                                 diffusionGeomCoeff=diffusionGeomCoeff,
     34                                                                 buildExplicitIfOther=buildExplicitIfOther)
     36     matrix += tmpMatrix
     37     RHSvector += tmpRHSvector

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\unaryTerm.py:62, in _UnaryTerm._buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
     46 """Build matrices of constituent Terms and collect them
     47 
     48 Only called at top-level by `_prepareLinearSystem()`
   (...)
     58 
     59 """
     61 if var is self.var or self.var is None:
---> 62     var, matrix, RHSvector = self._buildMatrix(var,
     63                                                SparseMatrix,
     64                                                boundaryConditions=boundaryConditions,
     65                                                dt=dt,
     66                                                transientGeomCoeff=transientGeomCoeff,
     67                                                diffusionGeomCoeff=diffusionGeomCoeff)
     68 elif buildExplicitIfOther:
     69     _, matrix, RHSvector = self._buildMatrix(self.var,
     70                                              SparseMatrix,
     71                                              boundaryConditions=boundaryConditions,
     72                                              dt=dt,
     73                                              transientGeomCoeff=transientGeomCoeff,
     74                                              diffusionGeomCoeff=diffusionGeomCoeff)

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\cellTerm.py:126, in CellTerm._buildMatrix(self, var, SparseMatrix, boundaryConditions, dt, transientGeomCoeff, diffusionGeomCoeff)
    123 b = numerix.zeros(var.shape, 'd').ravel()
    124 L = SparseMatrix(mesh=var.mesh)
--> 126 coeffVectors = self._getCoeffVectors_(var=var, transientGeomCoeff=transientGeomCoeff, diffusionGeomCoeff=diffusionGeomCoeff)
    128 dt = self._checkDt(dt)
    130 if inline.doInline and var.rank == 0:

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\cellTerm.py:86, in CellTerm._getCoeffVectors_(self, var, transientGeomCoeff, diffusionGeomCoeff)
     84 if self.coeffVectors is None or var is not self._var:
     85     self._var = var
---> 86     self._calcCoeffVectors_(var=var, transientGeomCoeff=transientGeomCoeff, diffusionGeomCoeff=diffusionGeomCoeff)
     88 return self.coeffVectors

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\cellTerm.py:69, in CellTerm._calcCoeffVectors_(self, var, transientGeomCoeff, diffusionGeomCoeff)
     68 def _calcCoeffVectors_(self, var, transientGeomCoeff=None, diffusionGeomCoeff=None):
---> 69     coeff = self._getGeomCoeff(var)
     70     weight = self._getWeight(var, transientGeomCoeff, diffusionGeomCoeff)
     71     if hasattr(coeff, 'old'):

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\term.py:480, in Term._getGeomCoeff(self, var)
    478 def _getGeomCoeff(self, var):
    479     if self.geomCoeff is None:
--> 480         self.geomCoeff = self._calcGeomCoeff(var)
    481         if self.geomCoeff is not None:
    482             self.geomCoeff.dontCacheMe()

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\sourceTerm.py:23, in SourceTerm._calcGeomCoeff(self, var)
     22 def _calcGeomCoeff(self, var):
---> 23     self._checkCoeff(var)
     25     if self.coeff.shape != () and self.coeff.shape[-1] != len(var.mesh.cellVolumes):
     26         return self.coeff[..., numerix.newaxis] * CellVariable(mesh=var.mesh, value=var.mesh.cellVolumes)

File ~\AppData\Roaming\Python\Python310\site-packages\fipy\terms\cellTerm.py:55, in CellTerm._checkCoeff(self, var)
     53 if var.rank == 0:
     54     if shape != ():
---> 55         raise TypeError("The coefficient must be rank 0 for a rank 0 solution variable.")
     57 if shape != () and len(shape) != 2 and shape[0] != shape[1]:
     58     raise TypeError("The coefficient must be a rank-0 or rank-2 vector or a scalar value.")

TypeError: The coefficient must be rank 0 for a rank 0 solution variable.

[guyer]

Your source S must be defined on the cells of the mesh (it must be a CellVariable).

Defining a source

S=(S_a/delta)*F*np.exp((-z/delta)-((4*math.log(2))*(((x*x)+(y*y))/(w0*w0))))*(1/tp)*math.sqrt((4*math.log(2)/math.pi))*np.exp((-4*math.log(2))*(((tt-2*tp)*(tt-(2*tp))/(tp*tp))))

in terms of 3D coordinates

x=np.linspace(0.1e-11,500e-9,100)
y=np.linspace(0.1e-11,500e-9,100)
z=np.linspace(0.1e-11,500e-9,100)

doesn't make any sense when you're using a 1D mesh.

Define your source in terms of the coordinates of m.

[Sandip433]

You can try this (solved in 1D mesh)

from fipy import *
import matplotlib.pyplot as plt
Lx=1e-6
dx=1e-9 
nx=Lx/dx 
x_0=Lx/2
T_e_t0=300.0 
G=1e17  
S_a=0.49
delta=15e-9
F=170
t_p=100e-15 
w_0=650e-9 
gamma=96.8  
C_l=3.62e6 
timeStepDuration=0.5e-15
K_e=391
K_l= 7
m=Grid1D(nx=nx, dx=dx)
x = m.cellCenters[0]
time = Variable(0.)
S=((S_a*F)/(t_p*delta)*numerix.sqrt((4*numerix.log(2)/numerix.pi))*numerix.exp(-x/delta)*numerix.exp(-4*numerix.log(2)*((x-x_0)/w_0)**2))*numerix.exp(-4*numerix.log(2)*((time-2*t_p)/t_p)**2)   
T_e=CellVariable(name="Electron Temperature", mesh=m, value=T_e_t0, hasOld=True)
T_l=CellVariable(name="Lattice Temperature", mesh=m, value=T_e_t0, hasOld=True)
eq_e=TransientTerm(coeff=gamma*T_e, var=T_e)==DiffusionTerm(coeff=K_e, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_e)+ImplicitSourceTerm(coeff=G, var=T_l)+S
eq_l=TransientTerm(coeff=C_l, var=T_l)==DiffusionTerm(coeff=K_l, var=T_l)+ImplicitSourceTerm(coeff=G, var=T_e)-ImplicitSourceTerm(coeff=G, var=T_l)
eq=eq_e & eq_l
results_T_e = []
results_T_l = []
t_list=[]
while time()<4e-12:
    t_list.append(float(time()))
    results_T_l.append(float(T_e.max()))
    results_T_e.append(float(T_l.max()))
    T_e.updateOld()
    T_l.updateOld()
    for i in range(1):
        res=eq.sweep(dt=timeStepDuration)
    time.setValue(time() + timeStepDuration)
plt.plot(t_list,results_T_e,'r',label="Electron Temperature")
plt.plot(t_list,results_T_l,'b',label="Lattice Temperature")
plt.xlabel('Time')
plt.ylabel('Temperature')
plt.legend(loc="upper right")
plt.show()