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: [email protected] <[email protected]> on behalf of Guyer,
Jonathan E. Dr. (Fed) <[email protected]>
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 <[email protected]> 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
[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 ]
