I have recently been trying to solve these 6 coupled, nonlinear PDEs of the
general form:
[cid:a9dfc2f3-d8a7-47c2-883b-b8466ce6d924]
(where N_1,N_2,P_1,P_2,E_1, and E_2 are the solutions sought; t and z are the
independent variables (only 2 at the moment); a,b,c,d,e,f,g,h and k are
positive constants)
I have set the parameters and constants equal to 1 for simplicity in this
example. Here is my code:
[code]
from fipy import *
import numpy as np
zmax = 10. #final length of z domain
tmax = 5. #final length of t domain
if __name__ == "__main__":
steps = nz = nt = 100 #number of elements, where dn is increments, so range
is nn*dn
#for the nth variable assessed
else:
steps = nz = nt = 10
mesh = Grid2D(nx=nz, ny=nt, dx=(zmax/nz), dy=(tmax/nt)) #dx = dz, dy = dt
z, t = mesh.cellCenters #pairs z and t to make a grid with nz by nt elements
dz = (zmax/nz)
dt = (tmax/nt)
N1 = CellVariable(name=r"$N_1$", mesh=mesh, hasOld=True, value=0.0) #value here
changes every element
N2 = CellVariable(name=r"$N_2$", mesh=mesh, hasOld=True) #and sets values of
function.old argument!
P1 = CellVariable(name=r"$P_1$", mesh=mesh, hasOld=True)
P2 = CellVariable(name=r"$P_2$", mesh=mesh, hasOld=True)
E1 = CellVariable(name=r"$E_1$", mesh=mesh, hasOld=True, value = 0.)
E2 = CellVariable(name=r"$E_2$", mesh=mesh, hasOld=True)
#The CellVariables are in a grid with nz by nt elements, which corresponds to
#a particular point on the grid
#The CellVariables are in the format of F = F(z,t) where t is constant for nz
elements
#and then changes for the nz + 1 element to the next step t + dt
for i in range(nz): #this sets initial condition at t = dt (not necessarily t =
0)
N1[i] = 0.01
N2[i] = 0.
P1[i] = 1.
P2[i] = 0.
E1[i] = 0.
E2[i] = 1.
#sets boundary condition at z = dz (not necessarily z = 0)
N1[::nt] = 0.01*(np.e**(-2.*np.log(2)*((t[::nt]-t[0])**2)))
N2[::nt] = 0.
P1[::nt] = 1.*(np.e**(-2.*np.log(2)*((t[::nt]-t[0])**2)))
P2[::nt] = 0.
E1[::nt] = 0.
E2[::nt] = 1.*(np.e**(-2.*np.log(2)*((t[::nt]-t[0])**2)))
if __name__ == "__main__":
viewer = Viewer(vars=(N1, N2, P1, P2, E1, E2)) # , datamin=0., datamax=1.)
#constants & parameters
hbar = 1
c = 1
eps = 1
omega = 1.
d = 1.
T1 = 1.
T2 = 1.
Ldiff = 1.
CoEff = CellVariable(mesh=mesh, value=(1), rank=1) #vector with values 1
dN1dt = (-2/hbar)*(P2*E1 + P1*E2) - N1/T1
dN2dt = -N2/T1
dP1dt = (2*(d**2)/hbar)*(E2*N1 - E1*N2) - P1/T2
dP2dt = (2*(d**2)/hbar)*(E1*N1 + E2*N2) - P2/T2
dE1dz = (omega/(2*eps*c))*P2 - E1/Ldiff
dE2dz = (omega/(2*eps*c))*P1 - E2/Ldiff
eq1 = (TransientTerm(var=N1) == dN1dt)
eq2 = (TransientTerm(var=N2) == dN2dt)
eq3 = (TransientTerm(var=P1) == dP1dt)
eq4 = (TransientTerm(var=P2) == dP2dt)
eq5 = (CentralDifferenceConvectionTerm(coeff=CoEff,var=E1) == dE1dz)
eq6 = (CentralDifferenceConvectionTerm(coeff=CoEff,var=E2) == dE2dz)
eq = eq1 & eq2 & eq3 & eq4 & eq5 & eq6
elapsedt = 0.
i = 0
while elapsedt < tmax:
N1.updateOld()
N2.updateOld()
P1.updateOld()
P2.updateOld()
E1.updateOld()
E2.updateOld()
i = i + 1
T = i*(tmax/nt)#min(100, numerix.exp(dexp))
elapsedt += dt
eq.solve(dt=T)
if __name__ == "__main__":
viewer.plot()
if __name__ == '__main__':
raw_input("Coupled equations. Press <return> to proceed...")
[\code]
I do not have any particular boundary/initial conditions I want to specifically
test at the moment, but besides trivial conditions, I have been finding the
output essentially always diverges. Increasing the number of steps also appears
to freeze the script. So please feel free to test any conditions you like even
if they differ from the ones in my code. Any thoughts on improving the code and
suggestions to solve the above nonlinear, coupled PDEs would be of much
interest. I am also curious if FiPy is well-suited to handle such a problem
involving nonlinear terms in the coupled PDEs since most examples I have noted
are linear and do not typically involve two independent variables.
Any and all input is welcome.
Abhilash
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
