Fair enough, thank you for the clarification. I've updated the code accordingly:
[code]
from fipy import *
import numpy as np
#constants & parameters
omega = 2.*np.pi*(1612.*10**(6)) #angular frequency of EM
eps = 8.85*10**(-12)
c = 3.*(10**8)
hbar = (6.626*10**(-34))/(2*np.pi)
eta = 0.01 #inversion factior
nn = 10.**7 #population density
n = eta*nn #inverted population density
lambdaOH = c/(1612.*10**(6)) #wavelength
gamma = 1.282*(10.**(-11))
Tsp = 1./gamma
TR = 604800.
T1 = 210.*TR
T2 = 210.*TR
L = (Tsp/TR)*(8.*np.pi)/((3.*(lambdaOH**2.))*n)
d = (3.*eps*c*hbar*(lambdaOH**2)*gamma/(4.*np.pi*omega))**0.5 #dipole
transition element
radius = np.sqrt(lambdaOH*L/np.pi)
A = np.pi*(radius**2)
V = (np.pi*(radius**2))*L
NN = n*V #total inverted population
constant3 = (omega*TR*NN*(d**2)/(2*c*hbar*eps*V))
Fn = 1. # Fresnel number = A/(lambdaOH*L)
Lp = 1.
Ldiff = Fn*L/0.35
theta0 = 4.7*(10**(-5))
zmax = L #final length of z domain
tmax = 500. #final length of t domain
if __name__ == "__main__":
steps = nz = nt = 500 #number of elements, where dn is increments, so range
is nn*dn
#for the nth variable assessed
else:
steps = nz = nt = 500
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 here changes
every element
P1 = CellVariable(name=r"$P_1$", mesh=mesh, hasOld=True) #and sets values of
function.old argument!
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)
N1.constrain(0.5*np.sin(theta0), where=mesh.facesBottom)
P1.constrain(0.5*np.cos(theta0), where=mesh.facesBottom)
P2.constrain(0., where=mesh.facesBottom)
E1.constrain(0., where=mesh.facesRight)
E2.constrain(0., where=mesh.facesRight)
ones = CellVariable(mesh=mesh, value=(1), rank=1) #vector with values 1
ones0 = CellVariable(mesh=mesh, value=(1), rank=0)
eq1 = (TransientTerm(var=N1) == ImplicitSourceTerm(coeff=-2.*E1, var=P2) +
ImplicitSourceTerm(coeff= -2.*E2, var=P1) + ImplicitSourceTerm(coeff=
-1./(T1/TR), var=N1))
eq2 = (TransientTerm(var=P1) == ImplicitSourceTerm(coeff=2.*E2, var=N1) +
ImplicitSourceTerm(coeff= -1./(T2/TR), var=P1))
eq3 = (TransientTerm(var=P1) == ImplicitSourceTerm(coeff=2.*E1, var=N1) +
ImplicitSourceTerm(coeff= -1./(T2/TR), var=P2))
eq4 = (CentralDifferenceConvectionTerm(coeff = ones, var=E1) ==
ImplicitSourceTerm(coeff=constant3, var=P2) +
ImplicitSourceTerm(coeff=-1./Ldiff, var=E1))
eq5 = (CentralDifferenceConvectionTerm(coeff = ones, var=E2) ==
ImplicitSourceTerm(coeff=constant3, var=P1) +
ImplicitSourceTerm(coeff=-1./Ldiff, var=E2))
eq = eq1 & eq2 & eq3 & eq4 & eq5
while res > 1e-10:
res = eq.sweep(dt=dt)
N1.updateOld()
P1.updateOld()
P2.updateOld()
E1.updateOld()
E2.updateOld()
print E1, E2
[/code]
When coupling the equations, should it be done separately for the partial
derivatives with respect to time and z (i.e. eq1, eq2, and eq3 are coupled
together, and eq4 and eq5 are coupled together since E1 or E2 is not updated as
N1, P1, and P2 are over the time steps)?
Also, with the current code, the variables do not appear to be evolving. I am
using a 2D mesh grid for both the temporal and spatial domain considered as I
would eventually like to see how N1, P1, P2, E1, and E2 vary both on z and t,
but is this the correct approach? It seems as though it might not be
appropriately handled by the CentralDifferenceConvectionTerm when doing so.
________________________________
From: fipy-boun...@nist.gov <fipy-boun...@nist.gov> on behalf of Guyer,
Jonathan E. Dr. (Fed) <jonathan.gu...@nist.gov>
Sent: June 27, 2016 3:59:30 PM
To: FIPY
Subject: Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
Strictly speaking, yes, you've couple the equations when you write
eq = eq1 & eq2 & eq3 & eq4 & eq5 & eq6
but you've not really couple anything because FiPy still only knows about the
diagonal elements of the coefficient matrix, so effectively, each equation is
being solved only for its "own" variable and not any of the others.
As far as FiPy is concerned,
TransientTerm(var=N1) = -2.*(P2*E1 + P1*E2) - N1/(T1/TR)
is an equation for N1, alone. FiPy does not see any implicit dependencies on
the right-hand side of the equation. Instead, do something like:
TransientTerm(var=N1) = ImplicitSourceTerm(coeff=-2. * E1, var=P2) +
ImplicitSourceTerm(coeff=-2. * E2, var=P1) + ImplicitSourceTerm(coeff=-1. /
(T1/TR), var=N1)
TransientTerm(var=P1) = ImplicitSourceTerm(coeff=2. * E2, var=N1) +
ImplicitSourceTerm(coeff=-1. / (T2/TR), var=P1)
TransientTerm(var=P2) = 2.*E1*N1 - P2/(T2/TR)
TransientTerm(var=P1) = ImplicitSourceTerm(coeff=2. * E1, var=N1) +
ImplicitSourceTerm(coeff=-1. / (T2/TR), var=P2)
This is only for the first three equations. You'll want to couple in E1 and E2
similarly.
Even when this is done, you'll still need to sweep, as your coefficients depend
on the other solution variables.
> On Jun 27, 2016, at 2:51 PM, Abhilash Mathews <amath...@uwo.ca> wrote:
>
> Thank you for the response. After looking at the second link specifically, is
> not the coupling formally performed in my inserted code by "eq = eq1 & eq2 &
> eq3 & eq4 & eq5 & eq6"? I made some changes to my code to be able to assess
> one case where I do know the initial/boundary conditions and have an expected
> output. This allowed me to make some simplifications and remove one equation
> altogether leading to 5 coupled equations:
>
>
>
> The conditions are:
>
> N1(z,0) = 0.5*np.sin(theta0)
> P1(z,0) = 0.5*np.cos(theta0)
> P2(z,0)= 0.
> E1(L,t) = 0.
> E2(L,t) = 0.
>
> where t = [0,500] and z = [0,L] where L is defined in the code.
>
> [code]
> from fipy import *
> import numpy as np
> import matplotlib.pyplot as plt
>
> #constants & parameters
> omega = 2.*np.pi*(1612.*10**(6))
> eps = 8.85*10**(-12)
> c = 3.*(10**8)
> hbar = (6.626*10**(-34))/(2*np.pi)
> eta = 0.01 #inversion factior
> nn = 10.**7 #population density
> n = eta*nn #inverted population density
>
> lambdaOH = c/(1612.*10**(6)) #wavelength
> gamma = 1.282*(10.**(-11))
> Tsp = 1./gamma
> TR = 604800.
> T1 = 210.*TR
> T2 = 210.*TR
> L = (Tsp/TR)*(8.*np.pi)/((3.*(lambdaOH**2.))*n)
> d = (3.*eps*c*hbar*(lambdaOH**2)*gamma/(4.*np.pi*omega))**0.5
> radius = np.sqrt(lambdaOH*L/np.pi)
> A = np.pi*(radius**2)
> V = (np.pi*(radius**2))*L
> NN = n*V #total inverted population
> constant3 = (omega*TR*NN*(d**2)/(2*c*hbar*eps*V))
> Fn = 1. # Fresnel number = A/(lambdaOH*L)
> Lp = 1.
> Ldiff = Fn*L/0.35
> theta0 = 4.7*(10**(-5))
>
> zmax = L #final length of z domain
> tmax = 500. #final length of t domain
> if __name__ == "__main__":
> steps = nz = nt = 50
> 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
> dz = (zmax/nz)
> dt = (tmax/nt)
>
> N1 = CellVariable(name=r"$N_1$", mesh=mesh, hasOld=True, value=0.0) #value
> here changes every element
> P1 = CellVariable(name=r"$P_1$", mesh=mesh, hasOld=True) #and sets values of
> function.old argument!
> 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)
>
> N1.constrain(0.5*np.sin(theta0), where=mesh.facesBottom)
> P1.constrain(0.5*np.cos(theta0), where=mesh.facesBottom)
> P2.constrain(0., where=mesh.facesBottom)
> E1.constrain(0., where=mesh.facesRight)
> E2.constrain(0., where=mesh.facesRight)
>
> if __name__ == "__main__":
> viewer = Viewer(vars=(N1, P1, P2, E1, E2)) # , datamin=0., datamax=1.)
>
> CoEff = CellVariable(mesh=mesh, value=(1), rank=1) #vector with values 1
>
> dN1dt = -2.*(P2*E1 + P1*E2) - N1/(T1/TR)
> dP1dt = 2.*E2*N1 - P1/(T2/TR)
> dP2dt = 2.*E1*N1 - P2/(T2/TR)
> dE1dz = constant3*P2 - E1/Ldiff
> dE2dz = constant3*P1 - E2/Ldiff
>
> eq1 = (TransientTerm(var=N1) == dN1dt)
> eq2 = (TransientTerm(var=P1) == dP1dt)
> eq3 = (TransientTerm(var=P2) == dP2dt)
> eq4 = (CentralDifferenceConvectionTerm(coeff=CoEff,var=E1) == dE1dz)
> eq5 = (CentralDifferenceConvectionTerm(coeff=CoEff,var=E2) == dE2dz)
>
> eq = eq1 & eq2 & eq3 & eq4 &
> elapsedt = 0.
> i = 0
>
> while elapsedt < tmax:
> N1.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)
>
> I_SR_scaled = (0.5*c*eps*(E1**2 + E2**2))
> I_SR = I_SR_scaled/((d*TR/hbar)**2)
> I_nc = (2./3.)*hbar*omega/(A*TR) #(should be) equivalent to I_nc
> IntensityRatio = (I_SR/(NN*I_nc))/0.001
>
> t = np.linspace(0,tmax,nt)
>
> IntensityRatioZequalL = []
> i = 0
> while i < (nt):
> A = IntensityRatio[(nz - 1) + i*nt]
> IntensityRatioZequalL.append(A)
> i = i + 1
>
> Intensity = np.array(IntensityRatioZequalL)
>
> fig = plt.figure(1)
> plt.plot(t,Intensity,'k:',linewidth = 1.5)
> plt.xlabel('t (ns)')
> plt.ylabel(r"$ \frac {I_{SR}}{NI_{nc}}$")
> plt.show()
>
> [/code]
>
> On top of the existing documentation, I've been following the example of
> cahnHiliard.mesh2Dcoupled closely and notice their mesh is x and y while I am
> assigning mine as z (I am only considering one spatial variable at the
> moment) and t. Without making the grid 2D, I was unable to set boundary
> conditions for z and t. (I assume the partial derivative with respect to z is
> set by the CentralDifferenceConvectionTerm term.) Although when running the
> code, it appears the z variable is not being integrated at each new step. I
> am also applying boundary conditions at the end of the sample interval (i.e.
> at z = L) which may affect how the variables are updated.
>
> Any thoughts on how to ensure both z and t are being updated appropriately in
> this case?
> From: fipy-boun...@nist.gov <fipy-boun...@nist.gov> on behalf of Guyer,
> Jonathan E. Dr. (Fed) <jonathan.gu...@nist.gov>
> Sent: June 27, 2016 10:45:55 AM
> To: FIPY
> Subject: Re: FiPy for nonlinear, coupled PDEs in multiple independent
> variables
>
> You have nothing here to actually couple the equations. At a minimum, you
> need to sweep each timestep to convergence. See:
>
>
> http://www.ctcms.nist.gov/fipy/documentation/FAQ.html#iterations-timesteps-and-sweeps-oh-my
>
> Better would be to formally couple the equations. See:
>
> http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#coupledequations
>
>
> > On Jun 16, 2016, at 2:48 PM, Abhilash Mathews <amath...@uwo.ca> wrote:
> >
> > I have recently been trying to solve these 6 coupled, nonlinear PDEs of the
> > general form:
> > <Screen Shot 2016-06-16 at 2.23.03 PM.png>
> > (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
> > fipy@nist.gov
> > http://www.ctcms.nist.gov/fipy
> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> <Screen Shot 2016-06-27 at 2.47.37
> PM.png>_______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
