Jonathan -
To estimate back diffusion and an accumulation at the interface (left-most
boundary), don't I want the precursor of the exponential to change at every
timestep as phi changes?
Forgive me if this should be evident, but want to make sure I understand.
Again, thank you for your assistance!!
________________________________
From: Jonathan Guyer <[email protected]>
To: FIPY list <[email protected]>
Sent: Monday, May 13, 2013 7:10 AM
Subject: Re: Transient Diffusion with Time Varying BC
On May 12, 2013, at 10:23 PM, Chuck Holbert <[email protected]> wrote:
> I'm looking for a bit of assistance in how to organize a FiPy script I am
> trying to write. I am trying to model transient 1D diffusion into and out of
> soil with an exponentially decaying boundary condition. I assuming a
> starting concentration and that diffusion occurs for an initial using a low
> first-order rate constant (-0.05/yr). After this period, an increased rate
> is applied (-1.0/yr) for three years and then the rate decreases back to the
> initial value (-0.05/yr).
>
> The script I have is below. I'm not sure the way I am setting my left
> boundary condition is correct. Any advice we be greatly appreciated.
I see two issues with the way you calculate the boundary condition:
valueLeft = phi[0] * numerix.exp(-1.0*(time-43))
The most significant one is that `phi[0]` returns a Variable quantity, so the
precursor of the exponential will change at every timestep as phi changes. You
can fix that by using `phi[0].value`.
The other issue is more a matter of generality. `phi[0]` happens to work
because this is a 1D problem and the Grid1D topology has the left-most cell in
index 0. FiPy is intended to deal with multidimensional domains and
unstructured meshes, so if you ever wanted to do this problem on a different
mesh, you'd find that `phi[0]` wouldn't work at all. I recommend always writing
it the "correct" way, even in 1D:
valueLeft = phi.faceValue[mesh.facesLeft.value].value *
numerix.exp(-1.0*(time-43))
i.e., you want the left-most face value and you want the static value at the
beginning of the time period. It would probably be slightly more accurate to
get the value from the previous boundary condition:
phi0 = phi.faceValue[mesh.facesLeft.value].value
del phi.faceConstraints
valueLeft = phi0 * numerix.exp(-1.0*(time-43))
If you calculate the prefactor after you delete the face constraints, then you
get the leftmost cell value, which is not quite the same as the leftmost face
value.
_______________________________________________
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 ]
