Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
On Tue, Jun 28, 2016 at 12:53 PM, Abhilash Mathews wrote: > > 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)) These are not great equations from FiPy's perspective as they are purely convective. However, if you want to persist with FiPy, you may want to use some sort of relaxation for the last two equations. One thing you could do is add in a TransientTerm and DiffusionTerm with small coefficients. Initially, just use coefficients of 1 to see if actually having those terms helps the numerical stability. If that helps, see if you can dial down those coefficients to get a more physical solution or use other techniques that maintain the correct physics. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
Thank you for the suggestions. I've updated the code accordingly with the same initial/boundary conditions previously mentioned and the graphing output will be inserted afterwards: [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. #not scaling z yet Ldiff = Fn*L/0.35 theta0 = 4.7*(10.**(-5.)) #initial Bloch angle = 2/np.sqrt(NN*eta) zmax = L/Lp #final length of z domain tmax = 500. #final length of t domain if __name__ == "__main__": steps = nz = 500 else: steps = nz = 50 mesh = Grid1D(nx=nz, dx=(zmax/nz)) z = mesh.cellCenters[0] dz = (zmax/nz) #where nz = steps in this case dt = tmax/steps N1 = CellVariable(name=r"$N_1$", mesh=mesh, value = 0.5*np.sin(theta0), hasOld=True) #value here changes every element P1 = CellVariable(name=r"$P_1$", mesh=mesh, value = 0.5*np.cos(theta0), hasOld=True) #and sets values of function.old argument! P2 = CellVariable(name=r"$P_2$", mesh=mesh, value = 0., hasOld=True) # N1(z,0) = 0.5sin(theta_0), P1(z,0) = 0.5cos(theta_0) E1 = CellVariable(name=r"$E_1$", mesh=mesh) #E1 and E2 are not transient terms E2 = CellVariable(name=r"$E_2$", mesh=mesh) #therefore hasOld != True E1.setValue(0., where = z > z[nz-2]) # E1(L,t) = 0 E2.constrain(0., where = z > z[nz-2]) # E2(L,t) = 0 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 res = 1. elapsedTime = 0 while elapsedTime < tmax: N1.updateOld() P1.updateOld() P2.updateOld() while res > 1e-10: res = eq.sweep(dt=dt) print res print N1, P1, P2, E1, E2 elapsedTime += dt [/code] (note: the boundary conditions are at the end of the sample [i.e. z = L]) After making the recommend changes, the code appears to be updating as intended, but the results are still diverging. To work around this, I was looking through the different solvers (e.g. trilinosNonlinearSolver) and was wondering if there are any in particular you think would be of help for this nonlinear problem? I was also considering modifying the step size based on the rate of change of the solutions, but with diverging output, am uncertain on the best way to go about that. To firstly yield finite results, do you see any improvements in the above code? From: fipy-boun...@nist.gov on behalf of Guyer, Jonathan E. Dr. (Fed) Sent: June 28, 2016 11:08:15 AM To: FIPY Subject: Re: FiPy for nonlinear, coupled PDEs in multiple independent variables > On Jun 27, 2016, at 6:44 PM, Abhilash Mathews wrote: > > Fair enough, thank you for the clarification. I've updated the code > accordingly: > 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)? You should couple all of the equations together, if you can. eq4 and eq5 are quasistatic, but the values of P1, P2, E1, and E2 are used in eq1, eq2, and eq3, so you want everything updating implicitly together. > > Also, with the current code, the variables do not appear to be evolving. This code is mixing up timesteps and sweeps. > while res > 1e-10: > res
Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
> On Jun 27, 2016, at 6:44 PM, Abhilash Mathews wrote: > > Fair enough, thank you for the clarification. I've updated the code > accordingly: > 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)? You should couple all of the equations together, if you can. eq4 and eq5 are quasistatic, but the values of P1, P2, E1, and E2 are used in eq1, eq2, and eq3, so you want everything updating implicitly together. > > Also, with the current code, the variables do not appear to be evolving. This code is mixing up timesteps and sweeps. > while res > 1e-10: > res = eq.sweep(dt=dt) > N1.updateOld() > P1.updateOld() > P2.updateOld() > E1.updateOld() > E2.updateOld() > print E1, E2 Sweeping is about achieving convergence on the non-linear elements of your equations at a given timestep. Once converged, you can then advance to the next timestep (using .updateOld()). You need two nested loops to achieve this. See the example at the end of: http://www.ctcms.nist.gov/fipy/documentation/FAQ.html#iterations-timesteps-and-sweeps-oh-my There are no TransientTerms for E1 and E2, so they should not be declared with `hasOld=True` and you should not call .updateOld() on them. This shouldn't be harmful, but in my experience it is sometimes. > 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. FiPy's meshes are purely spatial. They would not do the right thing if one of the dimensions is time. You would need to build up a separate 2D array if you want to visualize a sequence of time steps as a single image. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
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 on behalf of Guyer, Jonathan E. Dr. (Fed) 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/T
Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
d()
> 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 on behalf of Guyer,
> Jonathan E. Dr. (Fed)
> 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 wrote:
> >
> > I have recently been trying to solve these 6 coupled, nonlinear PDEs of the
> > general form:
> >
> > (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]
Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
The equations which appear to be missing from the previous email have been
attached.
From: fipy-boun...@nist.gov on behalf of Guyer,
Jonathan E. Dr. (Fed)
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 wrote:
>
> I have recently been trying to solve these 6 coupled, nonlinear PDEs of the
> general form:
>
> (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 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 p
Re: FiPy for nonlinear, coupled PDEs in multiple independent variables
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: fipy-boun...@nist.gov on behalf of Guyer,
Jonathan E. Dr. (Fed)
Sent: June 27, 2016 10:45:55 AM
To: FIPY
Subject: Re: FiPy for nonlinear, coupled PDEs in multiple independent variable
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 wrote:
>
> I have recently been trying to solve these 6 coupled, nonlinear PDEs of the
> general form:
>
> (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 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
FiPy for nonlinear, coupled PDEs in multiple independent variables
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 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 ]
