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:
[X]
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: [email protected] <[email protected]> on behalf of Guyer,
Jonathan E. Dr. (Fed) <[email protected]>
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 <[email protected]> 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
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
