The equations which appear to be missing from the previous email have been attached.
________________________________ 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 ]
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
