headers
Recif | 9 Jun 01:28
Favicon

Mixed Boundary Conditions


Hi,
 I would like to use FiPy on a problem with mixed (Robin) boundary conditions.
I was not clear on the implementation, so as a test I tried to solve the 
following problem for which the steady-state solution is known:
\begin{align*}
\frac{\partial C}{\partial t} &=\frac{\partial^2 C}{\partial x^2} 
  - Pe\frac{\partial C}{\partial x} -Da C \qquad 0 < x <1\\
x=0: Pe &= -\frac{\partial C}{\partial x} + Pe C \\
x=1: \frac{\partial C}{\partial x} & = 0
\end{align*}
My program is listed below. I had to switch the sign for the left boundary
condition from what I expected, to avoid negative concentrations. However, 
the solution does not converge to the analytical solution, and the flux on the 
right boundary is non-zero at later times, despite the specification of zero. 
I would appreciate any assistance in fixing these errors.

Recif
__________________________________________________
nx = 50
L=1.0
dx = L/nx
from fipy.meshes.grid1D import Grid1D
mesh = Grid1D(nx = nx, dx = dx)
from fipy.variables.cellVariable import CellVariable
C = CellVariable(name="concentration",mesh=mesh,value=0.0,hasOld=1)
Pe=1.0
Da=1.0
fluxLeft = -Pe*(C.getFaceValue()-1.0)
fluxRight = 0.0
(Continue reading)

Permalink | Reply |
headers
Daniel Wheeler | 11 Jun 18:37
Gravatar

Re: Mixed Boundary Conditions


Hi Recif,

I have debugged your code and eventually got the "right" answer.
Unfortunately, it required a heinous boundary condition hack to get it
to work. The script below works with trunk  (revision 2572) not with
version 1.2 and agrees with the analytical solution that you provided.

The boundary condition problems are due to boundary conditions being
applied to both convection terms and diffusion terms resulting in some
duplication. To get around these problems in your script I did two
things. One was to apply a fixed value boundary condition to the right
hand side, which has a value equal to the interior cell value. This
forces a zero derivative on both terms, which is what you want. The
second thing was to combine the convection and diffusion term on the
left side and see that the flux in should be -Pe. Unfortunately this
flux gets applied to both the diffusion and convection terms so I have
to divide it by 2. I'm sorry. This is really funky.

There are a number of things that come out of this example

  * Use this example as a test script documenting the boundary condition issues
  * As an intermediate fix for boundary conditions, have the fixedFlux
boundary condition apply only to the equation not to both the
diffusion term and convection term. I think we can do this relatively
easily.

In fact that intermediate fix may deal with 99% of  boundary condition
issues and we may not need to a attach boundary conditions to
variables at all.
(Continue reading)

Permalink | Reply |
headers
Daniel Wheeler | 12 Jun 17:56
Gravatar

Re: Mixed Boundary Conditions

I have fixed the issue with the fixedFlux boundary condition being
applied multiple times in revision 2573. Now a fixedFlux boundary
condition will apply to the whole equation rather than term by term.
There are still issues with boundary conditions, but this will help in
a number of cases. The new code for the Robin problem is below. I also
added it as a doctest that shows up in the manual. I've attached the
pdf for this.

Cheers

from fipy import *
nx = 100
dx = 1.0 / nx

mesh = Grid1D(nx=nx, dx=dx)
C = CellVariable(mesh=mesh)

D = 2.0
P = 3.0

BCs = (FixedFlux(faces=mesh.getFacesLeft(), value=-P),
       FixedValue(faces=mesh.getFacesRight(), value=C.getFaceValue()))

eq = PowerLawConvectionTerm((P,)) == \
     DiffusionTerm() - ImplicitSourceTerm(D)

A = numerix.sqrt(P**2 + 4 * D)
x = mesh.getCellCenters()[0]
CAnalytical = CellVariable(mesh=mesh)
CAnalytical.setValue(2 * P * exp(P * x / 2) * ((P + A) * exp(A / 2 * (1 - x))
(Continue reading)

Permalink | Reply |
headers
Recif | 18 Jun 15:24
Favicon

Re: Mixed Boundary Conditions


Hi,
 Based on the code you posted, I used the following program to test the
solution. It agrees with the steady-state solution, but not the time-varying
solution I get using a 1-D finite-volume code I wrote and the analytical
solution (I am unable to attach my graph using this access. I had sent a
subscription request but did not receive a response, so I cannot post to the
mailing list directly). I am using v2590 of FiPy from the trunk.
 Looking over the code, I am not clear about a couple of issues: 
1. Does the FixedFlux BC set the value of the total flux (sum of the convective
and diffusive fluxes) at the boundary? 
2. At the right-hand side boundary, FixedValue is used with a value of
C.getFaceValue(). Does this correspond to the cell value being set to the face
value there? The physical BC is for the diffusive flux to be zero at x=1.

Thanks for your assistance.

Recif

#!/usr/bin/env python
from fipy import *

nx = 100
L=1.0
dx = L/nx

mesh = Grid1D(nx = nx, dx = dx)
C = CellVariable(name="concentration",mesh=mesh,value=0.0,hasOld=1)

Da=1.0
(Continue reading)

Permalink | Reply |
headers
Daniel Wheeler | 18 Jun 17:17
Gravatar

Re: Mixed Boundary Conditions


On Wed, Jun 18, 2008 at 9:24 AM, Recif <recif@...> wrote:
>
> Hi,
>  Based on the code you posted, I used the following program to test the
> solution. It agrees with the steady-state solution, but not the time-varying
> solution I get using a 1-D finite-volume code I wrote and the analytical
> solution.

Can you post the transient analytical solution and I'll take a look.
Is the solution converging at each time step? The sweep loop just does
5 sweeps. It may need to do more. Use the residual as the stopping
criterion instead.

> (I am unable to attach my graph using this access. I had sent a
> subscription request but did not receive a response, so I cannot post to the
> mailing list directly).

Weird. Not sure why that would be the case. Maybe Jon knows.

> I am using v2590 of FiPy from the trunk.
>  Looking over the code, I am not clear about a couple of issues:
> 1. Does the FixedFlux BC set the value of the total flux (sum of the convective
> and diffusive fluxes) at the boundary?

It now applies to the whole equation (the sum). Previously it applied to both
the convection and diffusion terms, which was really broken.

> 2. At the right-hand side boundary, FixedValue is used with a value of
> C.getFaceValue(). Does this correspond to the cell value being set to the face
(Continue reading)

Permalink | Reply |
headers
Jonathan Guyer | 18 Jun 17:44
Favicon

Re: Mixed Boundary Conditions


On Jun 18, 2008, at 11:17 AM, Daniel Wheeler wrote:

>
>> (I am unable to attach my graph using this access. I had sent a
>> subscription request but did not receive a response, so I cannot  
>> post to the
>> mailing list directly).
>
> Weird. Not sure why that would be the case. Maybe Jon knows.

I'm afraid I don't. I have not received notice of your request, but  
I've now added you manually.
Permalink | Reply |
headers
Recif | 19 Jun 14:57
Favicon

Re: Mixed Boundary Conditions


Hi,
I have listed below the values of -Pe*C + grad C at faces 1 and 2 and
the values of grad C at faces n-1 and n for Pe=5. It seems that the
left hand BC is satisfied (=-Pe) at face 2 but not at face 1 (where
grad C is always reported as 0). Face n always shows a value of 0 for
grad C, which is consistent with the specified BC, but from the
behavior at face 1, I am not sure if the value reported for grad C at
face n is an initial value and not the actual calculated value. To
satisfy the BC at face n-1, even for a time step of 0.001 and 50
sweeps (using the value of res as a condition could result in an
infinite loop), grad C is about 0.01, which could be problematic for
mass conservation for long simulations. Also I am not clear on how to
specify the diffusive flux at a boundary to a non-zero value.

                           -Pe*C + grad C          grad C
sweeps=20 dt=0.01  t=0.1: -4.3292 -4.9779       -0.1983  0.  
          dt=0.01  t=0.5: -4.8049 -4.9903       -0.01868 0.
sweeps=50 dt=0.01  t=0.1: -4.3293 -4.9777       -0.05600 0.  
          dt=0.01  t=0.5: -4.8043 -4.9909       -0.01028 0.
          dt=0.001 t=0.1: -4.3199 -4.9770       -0.04070 0.
          dt=0.001 t=0.5: -4.8054 -4.9902       -0.00920 0.

Thanks.
Recif
Permalink | Reply |

Gmane