headers
Tony S Yu | 23 Jun 18:36
Favicon

Diffusion equation with a moving point source


Hi,

I'm trying to solve the diffusion equation with a moving point source.  
In otherwords,

	\partial\phi/\partial t = \nabla^2 \phi + \delta(x_s(t))

To represent the delta function on a discrete grid, I need to  
initialize a zero array and set certain points on the grid (given by  
position x_s) to non-zero values. Unfortunately, x_s is a function of  
time, but indexing an array with a Variable object doesn't retain lazy  
evaluation. Here's a toy example:

#############
from fipy import *

mesh = Grid1D(nx=10)
phi = CellVariable(mesh=mesh, value=0)
delta = CellVariable(mesh=mesh, value=0)
time = Variable(value=0)

x_s = time
x, = mesh.getCellCenters()
delta.setValue(1, where=(x > x_s) * (x < x_s + 1))

eq = TransientTerm() == ImplicitDiffusionTerm(coeff=1.) + delta
BCs = (FixedValue(faces=mesh.getFacesLeft(), value=0.),
        FixedValue(faces=mesh.getFacesRight(), value=0.))
for t in range(0, 10):
(Continue reading)

Permalink | Reply |
headers
Jonathan Guyer | 23 Jun 20:19
Favicon

Re: Diffusion equation with a moving point source


On Jun 23, 2008, at 12:38 PM, Tony S Yu wrote:

> In this example, the loop doesn't change `delta` because  
> `delta.setValue` evaluates the `x_s` Variable such that it's no  
> longer affected by changes to `time`.

Correct. setValue() is a one-time assignment.

> Is it possible to index arrays with Variable objects and retain lazy  
> evaluation capabilities?

I would do it like this:

xVar = CellVariable(mesh=mesh, value=x)
delta = 1. * ((xVar > x_s) & (xVar < x_s + 1))

In a future release of FiPy, x will, itself, be a CellVariable, but  
for the moment it's just an array, hence the need to declare xVar.

I multiply by 1. because delta otherwise consists of booleans, which,  
at a minimum changes the baseline of the solution. My guess is that  
something along the way is interpreting the booleans as 1 and -1 (or  
-1 and 0) instead of 1 and 0. I'm not seeing where, but multiplying by  
1. resolves the issue.

> Alternatively, if I move `delta.setValue` into the loop, I get  
> something similar to the desired solution (i.e. I have a moving  
> source).
(Continue reading)

Permalink | Reply |
headers
Tony S Yu | 24 Jun 01:35
Favicon

Re: Diffusion equation with a moving point source


On Jun 23, 2008, at 2:19 PM, Jonathan Guyer wrote:
> I would do it like this:
>
> xVar = CellVariable(mesh=mesh, value=x)
> delta = 1. * ((xVar > x_s) & (xVar < x_s + 1))

Thanks. This is exactly what I was trying to do.

>> But I believe this gives me an incorrect solution (checked against  
>> analytical solution) since the time integration would assume that  
>> the source has a fixed position during each time interval `dt`.
>
> The time integration is first order and the way I declare delta  
> above should have the exact same effect as moving delta.setValue  
> into the loop.

Hmm, that makes sense. Somehow I think I should have realized that  
these were the same.

> Do you get better agreement if you reduce the timestep?

Yeah, smaller timesteps help;  the convergence is a little slower than  
I imagined, though. I may have some other errors in my code that's  
causing problems (my actual code isn't the same as the example I gave).

Thanks for your help!
-Tony
(Continue reading)

Permalink | Reply |
headers
asreeve | 3 Jul 21:17

Setting up a rank-1 cell variable


Daniel and Jonathan,

I'm having some difficulty setting up a rank 1 cell variable that needs to 
be updated during sweeps.

I have three rank 0 cell variables:

head (dependant variable)
Top and Bot (values are fixed)

I calculate a saturation state (still rank 0):
satState=(1.*(head>=Top)+(head<Top)*(head-Bot)/(Top-Bot))

I need to use this to calculate a directional property.

CellSat=CellVariable(mesh=mesh,rank=1,value=1)
CellSat1=CellSat*[satState,satState,satState]

This seams to work but takes quite a long time (~2 minutes for 55000 
cells). Is there a better way to do this? I've also tried:

CellSat=CellVariable(mesh=mesh,rank=1,value=[satState,satState,satState])

but this also takes a long time.

Andy

----------------------
Andrew Reeve
(Continue reading)

Permalink | Reply |
headers
Jonathan Guyer | 3 Jul 21:58
Favicon

Re: Setting up a rank-1 cell variable


On Jul 3, 2008, at 3:17 PM, asreeve wrote:

> CellSat=CellVariable(mesh=mesh,rank=1,value=1)
> CellSat1=CellSat*[satState,satState,satState]
>
> This seams to work but takes quite a long time (~2 minutes for 55000  
> cells). Is there a better way to do this? I've also tried:
>
> CellSat 
> =CellVariable(mesh=mesh,rank=1,value=[satState,satState,satState])
>
> but this also takes a long time.

I can't imagine why. Both are virtually instantaneous for me, with a  
38**3 mesh. The latter is certainly the preferred syntax, though. The  
former sets up a pointless binOp that just multiplies everything by 1.  
Moreover, there's no need for the list of satStates;  
CellSat1=CellSat*satState should work just as well.

For that matter, you could just do

CellSat=CellVariable(mesh=mesh,rank=1,value=satState)

I guess I need to see more code.

Why does it need to be rank-1 if all the values are the same?
Permalink | Reply |

Gmane