Hi
Thank you for the solution. This indeed looks like a nifty way to solve the
issue at hand.
Being a novice in FiPy, and Finite Volume methods in general, I am a bit fuzzy
about a couple of things.
1. The backward Euler method is used in formulating the equation,
"n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) " . But Backward Euler is a
first order method, isn't it ? I am a bit confused about the second order
accuracy statement.
2. In, the equation, "eq = (fp.TransientTerm() ==
fp.DiffusionTerm(diffCoeff) - \
fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff *
mesh.facesRight).divergence)) " ,
from what I understand the 'ImplicitCoeff' is zeroed out in all other
locations other than the right boundary by the boolean method 'mesh.facesRight'
, right ?
If this is being implemented, do we have the restriction to use fixed dx ?
Although dx appears in the Implicit term, it can just be set to the value of
'dx' of the last node (i.e. length of the last segment in a variable-mesh sized
1D geometry file).Am I missing something here ?
3. Also, I am a bit confused about the negative signs used in the equation
above.The implicitCoeff is given a negative sign, and the equation's
ImplicitSourceTerm is also given a negative sign. From what I undertand,
this is what we intend to do ?
* Assign the diffusion coefficient D throughout the domain, using the
faceVariable method.
* Manually set it to zero at the right boundary
* Then **add** the Diffusion achieved by the div.(D grad( phi)) by using
the ImplicitSourceTerm method, and ensure that this effect gets added only at
the right boundary.
If this is the workflow, then perhaps we just have to declare the
ImplicitSourceTerm to be positive, and use the addition sign instead of
subtraction sign for the implicitSourceTerm in the equation ?
4. Yes, the k needs to be negative for the analytical solution to be stable.
My apologies. I should have been more clear in my earlier posts to the list.
5. Taking one massive step to get to the answer. Here is a big question
(due to my poor knowledge in the subject)
I understand that the implicitness enables us to take these large time-steps
due to the affine and A-stable properties of Implicit Methods.But are we
required to do this ?My problem is that, I am solving a coupled system of
PDEs and my other equations, require me to step very very slowly, due to the
inherent stiffness of my system. The PDEs describe processes of varying
time-constants. So, does taking the timestep to be small inhibit us in
any way ?
Once again, thanks a lot for your solution. Based on the discusssions with Ray
and from this particular method from your posting to the list, I think I
shall be able to solve the problem in 1D. The issues you pointed out about
using this approach in 2D went right over my head. I have to think about this,
when I implement them perhaps in a week or so.
Thanks once again,
Krishna
From: fipy-boun...@nist.gov on behalf of Daniel Wheeler
Sent: Friday, June 10, 2016 6:17 PM
To: Multiple recipients of list
Subject: Re: casting implicit Boundary Conditions in FiPy
Hi Krishina and Ray,
Thanks for the interesting discussion. I'm not 100% sure about everything that
Krishina is asking for in the latter part of the discussion so I'm just going
to address the code that Ray has developed below (my code is below). I think
there is a way to handle the right hand side boundary condition both implicitly
and with second order accuracy using an ImplicitSourceTerm boundary condition
trick. Sort of similar to what is here,
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-fixed-flux-boundary-conditions.
Like I said, I think the code below is both second order accurate (for fixed
dx) and implicit. Extending this to 2D might raise a few issues. The grid
spacing needs to be in the coefficient so, obviously, dx needs to be fixed in
both x and y directions. Also, I'm not sure if I'm missing a factor of dx in
the source term in 1D as I'm using both a divergence and an ImplicitSourceTerm
so there is some question about volumes and face areas in 2D as well. I'm also
confused about the signs, I had to flip the sign in front of the source to make
it work. It seems to look right though in 1D and you can just take one massive
time step to get to the answer. This will only work if k is negative otherwise
it's unstable, right?
The way I came up with the source was the following
n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2)
and then solve for n.grad(phi) which gives
n.grad(phi) = k * phi_P / (1 - dx * k / 2)
where phi_P is the value of phi at the cell next to the boundary with the
implicit boundary condition. Then we fake the outbound flux to be the
expression on the right of the equation.
Just as a general note it would be great in FiPy i
