Hi Dan,
Thank you for the clarifications. Indeed they were quite helpful in getting a
more concrete understanding of the issues
"I'm confused by your question. Backward Euler refers to the time stepping
scheme, which is indeed first order. The spatial discretization scheme is
second order and using the boundary condition trick that I demonstrated is
second order and also implicit."
** yes. What I am conveying is that, the spatial equation phi_at_boundary =
phi_at_lastcell + (gradient/flux at boundary)*(distance between the boundary
and last cell-node), is analogous to future_value = present_value + (slope
at future value * delta_t), a Backward Euler step in the time-domain.
n.grad(phi) = k * phi_at_boundary = k * (phi_at_lastcell + (gradient/flux at
boundary)*(distance between the boundary and last cell-node)), might be a
first order approximation ? Again, I am unsure of this claim. At a very first
glance, I thought the derivation was similar/analgous to Backward Euler
Well, now getting into real issue/question. I am trying to extend this concept
on a pseudo-2D domain. Let me explain.
1. The diffusion process happens only along the y-axis, ie. we have
anisotropic diffusion.
2. The same implicit Neumann BC applies for each x co-ordinate, i.e. we
have d_phi/dy along top surface = - k * phi at the top surface
My code is shown below.
I think that your derivation also holds good in this pseudo-2D diffusion case.
With that sorted out,
1. I substituted the 2D equivalent constructions instead of their 1D
counterparts
2. For this demo case, I chose a 2x2 grid, with just nx=2 and ny=2 and
unit spacing in each direction
3. Defined the diffusion coefficient as a rank 2 faceVariable, and set its
value to a suitable tensor, such as to restrict movement in 2D direction.
4. Constrain the diffusion coefficient along the top boundary to be zero.
This boundary is where the special implicit Neumann BC applies
5. Calculate the implicitCoeff with the same formula/code (since the
derivation is the same)
6. Define the source term coefficient as (mesh.faceNormals *
implicitCoeff * mesh.facesTop).divergence)
Now, my problem is the following. Although the simulation seems to yield a
plausible-looking viewer animation, showing diffusion along y-axis, the
mesh.faceNormals * implicitCoeff * mesh.facesTop).divergence return a
'AddoverFaceVariable' object of a much reduced dimension. The mesh.faceNormals
* implicitCoeff * mesh.facesTop correctly returns the expected vectors of size
12 (for this 2x2 problem, with unit spacing, there are 12 faceCenters). It
correctly sets only the top face to the implicit coefficient. But the
moment, the divergence operator is called upon, the size of the returned
'AddoverFaceVariable' object drops to 4, with the last 2 values being set to
the ImplicitCoeff. Also, it just becomes a 1D vector, with first two elements
being zero. Shouldn't this ImplicitSourceTerm coefficient remain as a 2D
matrix ?
The code runs, but I am worried if I am correctly implementing the boundary
condition. Is there any (recommended/standard) way to validate my solution ? I
don't have the analytical solution.
import fipy as fp
nx = 2
ny = nx
dx = 1.
dy = dx
mesh = fp.Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
phi = fp.CellVariable(name = "Field Variable",mesh = mesh)
valueBottom = 1.
phi.constrain(valueBottom, mesh.facesBottom)
D = 1.0
diffCoeff = fp.FaceVariable(mesh=mesh, rank=2, value=[[[0,0],[0,D]]])
diffCoeff.constrain(0., mesh.facesTop)
k = -1.0 # We are modelling the BC (d_phi/dx = k*phi) at the top surface
implicitCoeff = D*k / (1. - k*dx/2.) # in order to model implicit BC
# print (mesh.faceNormals * implicitCoeff * mesh.facesTop).divergence
eq = (fp.TransientTerm() == fp.DiffusionTerm(diffCoeff)
+ fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff *
mesh.facesTop).divergence))
viewer = fp.Viewer(vars=phi, datamin=0., datamax=1.)
viewer.plot()
############Transient Simulation ################################
timeStep = 0.9 * dx**2 / (2 * D)
steps = 75
for step in range(steps):
eq.solve(var=phi, dt=timeStep)
viewer.plot()
Krishna
-----Original Message-----
From: [email protected] [mailto:[email protected]] On Behalf Of Daniel
Wheeler
Sent: 13 June 2016 03:17
To: Multiple recipients of list <[email protected]>
Subject: Re: casting implicit Boundary Conditions in FiPy
On Sun, Jun 12, 2016 at 8:14 AM, Gopalakrishnan, Krishnakumar
<[email protected]<mailto:[email protected]>>
wrote:
>
> 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.
I'm confused by your question. Backward Euler refers to the time stepping
scheme, which is indeed first order. The spatial discretization scheme is
second order and using the boundary condition trick that I demonstrated is
second order and also implicit.
> 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 ?
Exactly, mesh.facesRight is just evaluates to a boolean array. There is also
mesh.exteriorFaces which gets you all exterior faces and then you can also do
logical operations based on physical location to get more specific masks to
capture boundaries or internal regions.
> If this is being implemented, do we have the restriction to use fixed dx ?
I guess not if you're careful, but it could be a bit fiddly getting the correct
distance.
> 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 ?
No, but you need to know how to get the correct distance. There are arrays of
cell to face distances, but I'd have to think about how to get that working
right.
> 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 ?
I was just flipping signs to make it work so I'm not sure. Maybe it would be
better to have them both positive and it would make more sense that way.
> 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.
You got it.
> 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 ?
Yup, I think so.
> 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 ?
Oh no. You can still take as small a time step as necessary to capture the
dynamics at the time scale of interest. It's just that you're not limited by
the horrible dx**2 / D time step limit associated with explicit diffusion.
Typically it's best to resolve to some dx based time scale (rather than dx**2)
that captures an interface being tracked or some feature of interest moving
around.
> 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 ?
The only down side to small time steps is that their small :-)
> 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.
Good luck with it.
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]<mailto:[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 ]
