Dear Dr. Wheeler,
Thank you for your quick answer! I've just tried to run the script with a float
time variable but it's all the same.
Here is the full script (I do not use the fipy viewer but pylab, make sure you
have it installed!).
Best regards,
Fanny
-----
>> Full Program
# H Permeation
import numpy as n
from scipy.special import erf
import pylab as plt
from fipy import *
## Mesh
X_scale= n.linspace(0,2.0e+4,1000)# grid 0-> 20 um ; units[nm]
Y_scale= n.zeros((shape(X_scale))) #zero concentration at t=0
Nbx=len(X_scale)
deltax= 1.0/float(Nbx)
mesh=Grid1D(nx=Nbx, dx = deltax)
print "Mesh: OK."
## Data
D = 2.0e+4 # diffusion coeff in [nm2/s] = D(cm2/s)*1.E+14
Cs = 1.0 # H concentration at entry surface
deltatemps= 1.0 # time step [s]
print "dt=", deltatemps, "s"
duration=15*60.0 # H charging test duration [s]
# Problem definition
cH=CellVariable(name="H concentration", mesh=mesh, value= Y_scale)
eq= TransientTerm() == DiffusionTerm(coeff = D) # Fick
BCs = (FixedValue(faces = mesh.getFacesRight(), value=0),
FixedValue(faces=mesh.getFacesLeft(),value=Cs))
#... and solving
elapsed =0
while elapsed <= duration:
eq.solve(var = cH, boundaryConditions = BCs, dt = deltatemps)
elapsed=elapsed+deltatemps
# Comparison with the erfc solution
Sol=1-erf(n.array(X_scale)/(2*sqrt(D*duration)))
# Graph plotting
Y_scale=tuple(cH)
plt.plot(X_scale,Y_scale, marker='x')
plt.plot(X_scale, Sol, marker='o')
plt.show()
-----Message d'origine-----
De : [email protected] [mailto:[email protected]] De la part de Daniel Wheeler
Envoyé : jeudi 11 août 2011 16:28
À : Multiple recipients of list
Objet : Re: time issues with implicit method
On Thu, Aug 11, 2011 at 9:56 AM, JAMBON Fanny 213198
<[email protected]> wrote:
>
> First of all I'm trying to simulate pure fickian diffusion through the plate;
> but though I've read and read again the documentation, I still haven't
> understood how the "real time during which the equation must be solved" has
> to be given to the solver, or put in the program.
You seem to have figured it out in your script. Just passing a "dt"
argument to solve is all it takes.
> The program I've written (end of this e-mail), is similar to the first mesh1D
> example of the documentation, but when comparing the hydrogen profile to the
> erfc solution, great discrepancies arise.
> It's as though the solver time went faster than the analytical solution time,
> since the solver has already attained the steady state when the hydrogen
> penetration is just at its beginning according the analytical solution.
Can you post the entire example so it runs. It's missing imports and
the like. Thanks.
> So, my question is: how should the "real time" duration be given to the
> solver and/or the program?
It seems like it should be correct
> - is "duration" equal to the number of iterations calling the 'eq.solve(...)'
> times the 'dt' given to the solver, which is what I understood from the
> documentation (as in the example mesh1D: t=Steps*TimeStepDuration, at least
> for explicit method); and in that case, where is my mistake ?
You have the right idea. Make sure you initialize the time variables
with a float rather than an integer, but I don't think that should be
a problem
> - or is it something else, and in that case, could someone please put me on
> the right track?
Post the entire script so that it runs and I'll try and figure it out.
--
Daniel Wheeler
