Thanks alot. I modified the script by adding a diffusion term, however it seems
not to converge into a solution. Is this given the boundary conditions?
I am not sure why it crashes.
About the complex term, I found in the manual that one can use conjugate(). In
the case of giving i*u^2, would it look like this:
TransientTerm(coeff=1., var=phi) * conjugate(phi) ?
The script without the complex part that crashes is shown below.
#!/usr/bin/env python
from fipy import *
from fipy import numerix
nx = 50
dx = 1. / float(nx)
mesh = Grid1D(nx=nx,dx=dx)
X = mesh.cellCenters[0]
phi = CellVariable(mesh=mesh, name="solution")
phi.setValue(0.5-0.5*numerix.tanh(10*(X-0.5)))
vi = Viewer(vars=phi,datamin=0.0, datamax=1.0)
vi.plot()
raw_input("Initialization ...")
phi.constrain(1., mesh.facesLeft)
phi.constrain(0., mesh.facesRight)
phi_sq = CellVariable(mesh=mesh)
phi_sq.setValue( phi*phi )
# We now represent the equation, where u = phi, du/dt + du/dx + 2 d^2u/dx^2 +
u^2 = 0 in fipy commands. du/dx : is a convection term with a unit scalar
coefficient, i.e. <SpecificConvectionTerm>(coeff=(1.,), # var=u) , du/dt : is a
transient term that one can treat as before, i.e. TransientTerm(var=u). one can
add a second order derivative as ExplicitDiffusionTerm(coeff=D)
eq = TransientTerm(coeff=1., var=phi) + ExponentialConvectionTerm(coeff=(1.,),
var=phi) + ExplicitDiffusionTerm(coeff=(2.,), var=phi) + phi_sq == 0.0
dt = 0.01
steps = 100
for step in range(steps):
eq.sweep(dt=dt)
#
phi_sq.setValue( phi * phi )
#
vi.plot()
raw_input("Press <return> ...")
Sergio Manzetti
[ http://www.fjordforsk.no/logo_hr2.jpg ]
[ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ]
Midtun
6894 Vangsnes
Norge
Org.nr. 911 659 654
Tlf: +47 57695621
[ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory
] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ]
From: "Thibault Bridel-Bertomeu" <[email protected]>
To: "fipy" <[email protected]>
Sent: Thursday, May 18, 2017 3:02:54 PM
Subject: Re: Problems running a simple PDE
A file ? You can export the result to a file yes, although I did not
demonstrate it in my script.
This is not fipy though, it’s about using numpy to write an array of values
into a file !
As for the comparison with an analytic solution, I invite you to take a look at
the mesh1D diffusion example where they create a phiAnalytical. To plot it, you
can do
vi = Viewer(vars=(phi, phiAnalytical))
vi.plot()
which should plot at the same time the solution of the iterative process from
fipy and the analytical value.
Cordialement,
--
T. BRIDEL-BERTOMEU
2017-05-18 14:46 GMT+02:00 Sergio Manzetti < [
mailto:[email protected] | [email protected] ] > :
OK, I see now. Where is this file you mention accessed? It does not come out in
the working directory.
Sergio Manzetti
[ http://www.fjordforsk.no/logo_hr2.jpg ]
[ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ]
Midtun
6894 Vangsnes
Norge
Org.nr. 911 659 654
Tlf: [ tel:+47%2057%2069%2056%2021 | +47 57695621 ]
[ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory
] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ]
From: "Thibault Bridel-Bertomeu" < [ mailto:[email protected] |
[email protected] ] >
To: "fipy" < [ mailto:[email protected] | [email protected] ] >
Sent: Thursday, May 18, 2017 2:47:24 PM
Subject: Re: Problems running a simple PDE
The PDE you want to solve has a transient term, which is why the function
evolves in time.
However, as you can see it eventually converges to a given shape and stay in
that shape : this « stabilized shape » corresponds to the steady-state result
of your equation.
All you can see before are the transient organizations.
As I said in my previous e-mail, you can access phi as a numpy array and do
with it as you want : this phi is the discrete representation of your
steady-state result.
Cordialement,
--
T. BRIDEL-BERTOMEU
2017-05-18 14:40 GMT+02:00 Sergio Manzetti < [
mailto:[email protected] | [email protected] ] > :
BQ_BEGIN
I am just not sure on how to understand the actual beautiful result from the
PDE.
When initializing is pressed, it animates into a wave-like pattern.
I thought the numerical solution to the PDE was an actual function, but this
animated plot must be some form of time-dependent changing function? How can it
be given as a function in addition to the graph?
Sergio Manzetti
[ http://www.fjordforsk.no/logo_hr2.jpg ]
[ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ]
Midtun
6894 Vangsnes
Norge
Org.nr. 911 659 654
Tlf: [ tel:+47%2057%2069%2056%2021 | +47 57695621 ]
[ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory
] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ]
From: "Thibault Bridel-Bertomeu" < [ mailto:[email protected] |
[email protected] ] >
To: "fipy" < [ mailto:[email protected] | [email protected] ] >
Sent: Thursday, May 18, 2017 2:29:23 PM
Subject: Re: Problems running a simple PDE
I invite you to take a look at : [
http://www.ctcms.nist.gov/fipy/examples/README.html |
http://www.ctcms.nist.gov/fipy/examples/README.html ] for an overview of all
the examples written for fipy
and [ http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html |
http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html ] for a
reminder on how the finite volume method works.
As for your equation, if you are in 1D, you may consider it as follow :
du/dx : is a convection term with a unit scalar coefficient, i.e.
<SpecificConvectionTerm>(coeff=(1.,), var=u) - you will find a list of all
available convection schemes here : [
http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#convection-term
|
http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#convection-term
]
du/dt : is a transient term that you can treat as before, i.e.
TransientTerm(var=u)
and u^2, as a nonlinear term, one possibility is to update it during the
iterations.
All in all, you would write for instance :
#!/usr/bin/env python
from fipy import *
from fipy import numerix
nx = 50
dx = 1. / float(nx)
mesh = Grid1D(nx=nx,dx=dx)
X = mesh.cellCenters[0]
phi = CellVariable(mesh=mesh, name="solution")
phi.setValue(0.5-0.5*numerix.tanh(10*(X-0.5)))
vi = Viewer(vars=phi,datamin=0.0, datamax=1.0)
vi.plot()
raw_input("Initialization ...")
phi.constrain(1., mesh.facesLeft)
phi.constrain(0., mesh.facesRight)
phi_sq = CellVariable(mesh=mesh)
phi_sq.setValue( phi*phi )
eq = TransientTerm(coeff=1., var=phi) + ExponentialConvectionTerm(coeff=(1.,),
var=phi) + phi_sq == 0.0
dt = 0.01
steps = 100
for step in range(steps):
eq.sweep(dt=dt)
#
phi_sq.setValue( phi * phi )
#
vi.plot()
raw_input("Press <return> ...")
The sweep method allows to advance the equation of one timestep only. Then I
can update phi_sq which the next sweep will use to solve the equation. And so
on ..
Hope this helps to understand. I can however only advise you to browse through
all the examples, you will definitely find something similar to your case you
can start from !
Best
T. Bridel-Bertomeu
2017-05-18 14:10 GMT+02:00 Sergio Manzetti < [
mailto:[email protected] | [email protected] ] > :
BQ_BEGIN
OK, this really helped. If I wanted to change this into:
du/dx + i*du/dt + u^2 = 0
in the following format:
eqX = TransientTerm(var=phi) == ExplicitDiffusionTerm(coeff=D, var=phi)
for step in range(steps):
eqX.solve(dt=timeStepDuration)
there would be alot of different terms. Where can I find an explanation on how
to change these variables you mentioned into an equation one requires to run on
fipy?
Sergio
Sergio Manzetti
[ http://www.fjordforsk.no/logo_hr2.jpg ]
[ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ]
Midtun
6894 Vangsnes
Norge
Org.nr. 911 659 654
Tlf: [ tel:+47%2057%2069%2056%2021 | +47 57695621 ]
[ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory
] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ]
From: "Thibault Bridel-Bertomeu" < [ mailto:[email protected] |
[email protected] ] >
To: "fipy" < [ mailto:[email protected] | [email protected] ] >
Sent: Thursday, May 18, 2017 2:09:47 PM
Subject: Re: Problems running a simple PDE
Hello again Sergio,
When you define eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D)
you ready the equation, and then eqX.solve(var=phi, dt=timeStepDuration)
solves the equation for an infinitesimal timestep dt, and the variable phi as
asked.
You could also have written ::
eqX = TransientTerm(var=phi) == ExplicitDiffusionTerm(coeff=D, var=phi)
for step in range(steps):
eqX.solve(dt=timeStepDuration)
Best
Thibault Bridel-Bertomeu
2017-05-18 13:57 GMT+02:00 Sergio Manzetti < [
mailto:[email protected] | [email protected] ] > :
BQ_BEGIN
Dear all, referring to the email below, I removed the phi_nw and phi_old,
however, I am not sure where the actual PDE given in the first example at . 58,
d phi/dt = A d^2x∕dx^2
appears in the actual python code.
This would be important to change the PDE into other PDEs
Where do I define this PDE?
Thanks
Sergio Manzetti
[ http://www.fjordforsk.no/logo_hr2.jpg ]
[ http://www.fjordforsk.no/ | Fjordforsk AS ] [ http://www.fjordforsk.no/ | ]
Midtun
6894 Vangsnes
Norge
Org.nr. 911 659 654
Tlf: [ tel:+47%2057%2069%2056%2021 | +47 57695621 ]
[ http://www.oekolab.com/ | Økolab ] | [ http://www.nanofact.no/ | Nanofactory
] | [ http://www.aq-lab.no/ | AQ-Lab ] | [ http://www.phap.no/ | FAP ]
From: "sergio manzetti" < [ mailto:[email protected] |
[email protected] ] >
To: "fipy" < [ mailto:[email protected] | [email protected] ] >
Sent: Thursday, May 18, 2017 1:36:57 PM
Subject: Problems running a simple PDE
Hello, following the manual, at the first example with the 1D diffusion, I have
tried to make the python script for this (python2.7) and I get a missing name,
which is not defined in the example at all.
This is the name "cell", which appears on p 58 in the manual.
Please see this code which describes this example:
Thanks!
#Python script for solving a Partial Differential Equation in a diffusion case
from fipy import *
nx = 50
dx = 1
mesh = Grid1D(nx=nx, dx=dx)
phi = CellVariable(name="Solution variable",
mesh=mesh,
value=0.)
#We'll use the coefficient set to D=1
D=1
# We give then 2000 time-steps for the computation, which is defined by 90
percent of the maximum stable timestep, which is given
timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 2000
#Then we define the boundary conditions
valueLeft = 1
valueRight = 0
#The boundary conditions are represented as faces around the exterior of the
mesh, so we define the constraint by the given values on the left and right
side of the boundary using the phi.constrain() command
phi.constrain(valueLeft, mesh.facesRight)
phi.constrain(valueRight, mesh.facesLeft)
# We can also omit giving boundary conditions, then the default bc is
equivalent to a zero gradient.
#At this stage we define the partial differential equation as a numerical
problem to be solved
for step in range(steps):
for j in range(cells):
phi_new[j] = phi_old[j] \
+ (D * dt / dx**2) * (phi_old[j+1] - 2 * phi_old[j] + phi_old[j-1])
time += dt
#and additional code for the boundary conditions
eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D)
#We then want to view the results of the calculation and use the command
Viewer() for this
phiAnalytical = CellVariable(name="analytical value",
mesh=mesh)
if __name__ == '__main__':
viewer = Viewer(vars=(phi, phiAnalytical),
datamin=0., datamax=1.)
viewer.plot()
#If we have a semi-infinite domain, then the solution for this PDE is
phi=1-erf(x/2(Dt)^(0.5)). This requires to import SciPy library, so we import
that.
x = mesh.cellCenters[0]
t = timeStepDuration * steps
try:
from scipy.special import erf
phiAnalytical.setValue(1-erf(x / (2 * numerix.sqrt(D * t))))
except ImportError:
print ("The Scipy Library is not avaliable to test the solution to the given
trasient diffusion equation")
# The equation is then solved by repeatedly looping
for step in range(steps):
eqX.solve(var=phi,
dt=timeStepDuration)
if __name__ == '__main__':
viewer.plot()
print (phi.allclose(phiAnalytical, atol = 7e-4))
1
if __name__ == '__main__':
raw_input("Explicit transient diffusion. Press <return> to proceed...")
_______________________________________________
fipy mailing list
[ mailto:[email protected] | [email protected] ]
[ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ]
[ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy |
https://email.nist.gov/mailman/listinfo/fipy ] ]
_______________________________________________
fipy mailing list
[ mailto:[email protected] | [email protected] ]
[ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ]
[ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy |
https://email.nist.gov/mailman/listinfo/fipy ] ]
_______________________________________________
fipy mailing list
[ mailto:[email protected] | [email protected] ]
[ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ]
[ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy |
https://email.nist.gov/mailman/listinfo/fipy ] ]
BQ_END
_______________________________________________
fipy mailing list
[ mailto:[email protected] | [email protected] ]
[ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ]
[ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy |
https://email.nist.gov/mailman/listinfo/fipy ] ]
_______________________________________________
fipy mailing list
[ mailto:[email protected] | [email protected] ]
[ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ]
[ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy |
https://email.nist.gov/mailman/listinfo/fipy ] ]
BQ_END
_______________________________________________
fipy mailing list
[ mailto:[email protected] | [email protected] ]
[ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ]
[ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy |
https://email.nist.gov/mailman/listinfo/fipy ] ]
_______________________________________________
fipy mailing list
[ mailto:[email protected] | [email protected] ]
[ http://www.ctcms.nist.gov/fipy | http://www.ctcms.nist.gov/fipy ]
[ NIST internal ONLY: [ https://email.nist.gov/mailman/listinfo/fipy |
https://email.nist.gov/mailman/listinfo/fipy ] ]
BQ_END
_______________________________________________
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 ]
