Front-fixing method for moving boundary problem

Hi Jonathan and Daniel,

Jonathan, although I'm starting a new thread here, I first want to thank 
you for your once again very complete answer on the use of operator 
variables to get access to additional functions. I'm not 100\% sure yet 
whether I'll need to use one of those. Nevertheless, I found your answer 
extremely useful, it helped me a lot to understand more of the mechanics 
of FiPy.

Now, I'd like to have your thoughts on another aspect of my problem 
which is basically a Stefan problem or a free or moving boundary 
problem, whatever you like to call it. I want to use a front-fixing 
method to alleviate the moving boundary difficulty. Basically, I do a 
coordinate transformation to get a non-dimensional space coordinate, 
which allows me to use regular boundary conditions. However, this 
results in an additional ODE that needs to be solved along with the PDEs.

Here is an example, not the problem I want to solve, but the idea is the 
same:

\[
\frac{\partial^2 u}{\partial \xi ^2} = s^2 \frac{\partial u}{\partial t}
- s \xi \frac{ds}{dt}\frac{\partial u}{\partial \xi}
\]

where $\xi$ is the non-dimensional space coordinate, $s = s(t)$ and
\[
-\frac{1}{s} \left. \frac{\partial u}{\partial \xi}\right|_{\xi = 1}
= \frac{ds}{dt}

Re: Problems on Installing New Fipy

On Jul 16, 2008, at 5:34 PM, Edwin Garcia wrote:

> I think pylab is a module that FiPy uses for one of the viewers. At  
> least in the version that I have, it *is* part of the installation  
> and it is used to set the color bar. Mike can give you more details,  
> and even some comments on things that he has been modifying in his  
> installation to get the effects that he is trying to create.

I know what pylab is. We use it extensively. pylab supplies `ion()`,  
so the error message "'module' object has no attribute 'ion'" tells me  
that pylab is broken somehow. The most likely cause is that there is  
some file or module on Zhiwen's PYTHONPATH that is named "pylab" for  
some reason.

Problems on Installing New Fipy

faceVariable missing copy method

Hi,

I think `faceVariable` may be missing a `copy` method. I've included a  
patch below, which is basically a copy of `cellVariable`'s `copy`  
method (minus it's `hasOld` attribute and changing the rename from +  
'_old' to + '_copy').

Also, the patch includes a fix for matplotlib2DViewer.py, which has  
changed since I originally submitted the patch.

Cheers,
-Tony

In case you're curious about a use case: I wanted to set *almost* all  
of the exterior faces to a no flux boundary condition. So I naturally  
I tried

 >>> no_flux_faces = mesh.getExteriorFaces()
 >>> no_flux_faces.setValue(array[False], where=faces_with_other_BC)

But it turns out the above commands change the faces that the mesh  
recognizes as exterior (`getExteriorFaces` doesn't automatically copy;  
for performance reasons, maybe?). So I just tried to copy before  
modifying:

 >>> exterior = mesh.getExteriorFaces()
 >>> no_flux_faces = exterior.copy()

But, `faceVariable`'s `copy` method defaults to `Variable`'s `copy`  
method. As a result, the copied variable doesn't have a `mesh`  
attribute.

Re: Programing Differential Equations in Fipy

Did you install this version of fipy from your working copy?

Execute these commands at the python command line and tell me what happens.

    >>> import fipy, os
    >>> print fipy.__version__
    >>> print fipy.__path__
    >>> os.system('gmsh --version')

On Fri, Jul 11, 2008 at 4:13 PM, Zhiwen Liang <liangz@...> wrote:
> Hi Daniel,
> I downloaded the current version of Fipy using SVN. I just tried the
> diffusion example "circle.py" and it seems there is a problem with the GMSH
> mesh. The error is:
>   File "circle.py", line 301, in ?
>     exec(fipy.tests.doctestPlus._getScript())
>   File "<string>", line 39, in ?
>   File
> "/Users/liangz/stow/Fipy_current//lib/python2.3/site-packages/fipy/meshes/numMesh/gmshImport.py",
> line 394, in __init__
>     mesh2D.Mesh2D.__init__(self, **_DataGetter().getData(filename,
> dimensions = 2, coordDimensions = coordDimensions))
>   File
> "/Users/liangz/stow/Fipy_current//lib/python2.3/site-packages/fipy/meshes/numMesh/gmshImport.py",
> line 263, in getData
>     vertexCoords = self._calcVertexCoords(coordDimensions)
>   File
> "/Users/liangz/stow/Fipy_current//lib/python2.3/site-packages/fipy/meshes/numMesh/gmshImport.py",
> line 286, in _calcVertexCoords
>     numVertices = int(self.inFile.readline())
> ValueError: invalid literal for int(): 2 0 8
> You know what is the problem?
> Thanks,
> Zhiwen
> On Wed, Jul 9, 2008 at 1:19 PM, Daniel Wheeler <daniel.wheeler2@...>
> wrote:
>>
>> On Wed, Jul 9, 2008 at 11:18 AM, Zhiwen Liang <liangz@...> wrote:
>> > Hi Daniel,
>> > Thank you very much for providing such useful information! I really
>> > appreciated!
>> > My version is old and does not have the anisotropic diffusion
>> > coefficient. I
>> > will try the new version. This feature looks very powerful!
>>
>> Update to trunk at this point. We still haven't been able to do a release.
>>
>> > Regarding my second question, I did mean solving two equations with two
>> > unknowns. Is that possible?
>>
>> Yes, include both equations in an inner loop at each timestep and
>> iterate to convergence . Take a look at
>> <http://matforge.org/fipy/browser/trunk/examples/phase/binary.py#L636>
>> for an example.
>>
>>
>>
>> --
>> Daniel Wheeler
>
>

--

-- 
Daniel Wheeler

Re: Error in code

On Jul 10, 2008, at 9:58 AM, Rashmi Ranjan Mohanty wrote:

> I am running a Fipy code with a system size of 400x400 mesh points  
> on a Mac Pro 3GHz Dual-Intel Xeon processor with 2GB memory. After  
> running few steps it stops with the following error message:
>
> python(1013) malloc: *** vm_allocate(size=7921664) failed (error  
> code=3)
> python(1013) malloc: *** error: can't allocate region
> python(1013) malloc: *** set a breakpoint in szone_error to debug
> malloc fails for local dworkptr[].Traceback (most recent call last):
>  File "isothermal-new2D.py", line 355, in <module>
>    compRes = compEq.sweep(var = comp, boundaryConditions = BCsC,  
> dt=time, solver=solver)
>  File "/sw/lib/python2.5/site-packages/fipy/terms/term.py", line  
> 153, in sweep
>    self._solveLinearSystem(var, solver, matrix, RHSvector)
>  File "/sw/lib/python2.5/site-packages/fipy/terms/term.py", line  
> 102, in _solveLinearSystem
>    solver._solve(matrix, array, RHSvector)
>  File "/sw/lib/python2.5/site-packages/fipy/solvers/ 
> linearLUSolver.py", line 92, in _solve
>    LU = superlu.factorize(L._getMatrix().to_csr())
> MemoryError
>
> I guess it has something to do with memory, but could not figure out  
> a solution.
> I would appreciate your help.

My suspicion would lie on your use of the LU solver. The PySparse LU  
solver is robust, and often faster than the iterative solvers, but it  
does so by using a great deal of memory.

Have you tried other solvers?

Setting boundary conditions on interior cells

This is probably a bad idea, but I want to set a boundary condition on  
an interior cell. Alternatively, I would just have a rectangular mesh  
with some of the nodes masked off (kind of like a step on one of the  
rectangle's sides). The problem with changing the mesh in this manner  
is that this boundary condition needs to move in time.

My first thought was to force a cell to a given value by repeatedly  
resetting its value inside the solution loop. But, since the diffusion  
terms use the face values (not cell values) in the solution, I'm not  
getting the value back. The next (bad) idea I had was to set the face  
values, but I didn't see a way of doing this.

Is there a way of setting a boundary condition on an interior cell?

Thanks,
-Tony

Test script:
~~~~~~~~~~~~~~~
from fipy import *

dt = 0.001
# material properties
D = 1.
c0 = 1.
c_sink = 0.
# spacial parameters
nx = ny = 3
dx = dy = 1.
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
x, y = mesh.getCellCenters()
bottom_center = (x == 1.5) & (y == 0.5)

phi = CellVariable(name='concentration', mesh=mesh, value=c0)
phi.setValue(c_sink, where=bottom_center)
eq = TransientTerm() == ImplicitDiffusionTerm(coeff=D)

all_faces = mesh.getExteriorFaces()
BCs = (FixedFlux(faces=all_faces, value=0),)

for step in range(3):
     phi.setValue(c_sink, where=bottom_center)
     eq.solve(var=phi, boundaryConditions=BCs, dt=dt)
     # one of the cells should be 0 (c_sink value)
     print phi.reshape((ny, nx))

Re: Programing Differential Equations in Fipy

On Wed, Jul 9, 2008 at 11:18 AM, Zhiwen Liang <liangz@...> wrote:
> Hi Daniel,
> Thank you very much for providing such useful information! I really
> appreciated!
> My version is old and does not have the anisotropic diffusion coefficient. I
> will try the new version. This feature looks very powerful!

Update to trunk at this point. We still haven't been able to do a release.

> Regarding my second question, I did mean solving two equations with two
> unknowns. Is that possible?

Yes, include both equations in an inner loop at each timestep and
iterate to convergence . Take a look at
<http://matforge.org/fipy/browser/trunk/examples/phase/binary.py#L636>
for an example.

--

-- 
Daniel Wheeler

(unknown)

Daniel and Jonathan:

I'd like to concatenate two 3-D extruded mesh files. The examples in the 
users guide indicate this can be done, but I can't seem to get the 
examples listed there to work and was curious if this feature has been 
dropped or the convention for doing this has changed?

I've created two grids that I'd like to add together in the Z direction. 
Some of the cells have a thickness of zero, and need to be removed from 
the mesh. I was hoping I could do this be extruding a mesh that removed 
the 'empty cells', and then attach it to the other mesh. Is this possible, 
or do I need to handle all this in one extruded file?

Andy

----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469

Re: Programing Differential Equations in Fipy

Hi Zhiwen,

Hope all is well.

On Mon, Jul 7, 2008 at 4:44 PM, Zhiwen Liang <liangz@...> wrote:
> Hi all,
>
> Professor Garcia and I are working on solving some differential equations
> that look different from those that can be assembled by the regular Fipy
> terms. It would be great if you can give us some ideas.
>
> First, we are dealing with terms that are cross differentials with respect
> to different axis. For example, \frac{ \partial^2 \phi }{\partial x \partial
> y}.
>
> Second, the variables that are being solved for are embedded in several
> equations. For example, we want to solve these two equations: (Sorry they
> are not good examples.)
> \frac{\partial^2 \phi_1}{\partial x^2}+\frac{\partial^2 \phi_2}{\partial
> y^2}=0
> \frac{\partial^2 \phi_1}{\partial x \partial y}+\frac{\partial^2
> \phi_2}{\partial x \partial y}=0

We now have anisotropic diffusion, which, I think, deals with the
terms above. For example,

    \frac{\partial^2 \phi_1}{\partial x \partial y}

can be represented by,

    tensor_coeff = ((0, 1), (0, 0))
    DiffusionTerm((tensor_coeff,))

tensor_coeff refers to \Gamma_{ij} in the operator

          \partial_i  \Gamma_{ij} \partial_j

for example.

Basically, the diffusion coefficient can now be a tensor. This gives a
lot more flexibility in the equations that can be defined. We only
have one example of this usage thus far
<http://matforge.org/fipy/browser/trunk/examples/diffusion/anisotropy.py>.

For your equations above, it might be worth transforming the variables
to something like, psi_1 = phi_1 + phi_2 and psi_2 = phi_1 - phi_2.

BTW I'm assuming you are not asking a question about solving two
equations and two unknowns. Correct?

> And I wonder if there are people trying to solve for elasticity problems
> using Fipy.

There has been some chatter, but I am not sure how successful attempts
have been <http://search.gmane.org/?query=elasticity&group=gmane.comp.python.fipy>.

Cheers

--

-- 
Daniel Wheeler

Simulating complex PDEs

Hi,

A beginner's question: I have a set of coupled partial differential
equations which I feel that I am close to dealing with in FiPy,
except that they are in phasor form (that is, they have complex-valued
variables).  When I try introducing a source term with complex-typed
coefficient or a variable, initialized with a non-zero imaginary
part, I get this error:

   TypeError: array cannot be safely cast to required type

Is FiPy capable of solving such equations in their complex form?
If not, is the best alternative just splitting each equation into
a pair of coupled equations for the real and imaginary parts of the
complex variable?  I'm guessing that this would be essentially
equivalent to solving any other sets of coupled differential
equations.

(For concreteness, the equations are roughly of this form:

eq1 = TransientTerm(coeff = 1) - VanLeerConvectionTerm(coeff =
(-vg,),diffusionTerm = diffTerm) -  ImplicitSourceTerm(coeff =
-alpha) - ImplicitSourceTerm(coeff=gB*Ip)

except that the (constant) Ip in the last source term will have to
be changed to the product of two other variables.)

Any advice is appreciated.

Thanks!
David