Re: error in gradient even with orthogonal grid

[Daniel Wheeler]

Hi Jamie,

Sorry for the delayed reply.

On Thu, Jul 21, 2016 at 11:31 AM, James Pringle  wrote:
>
> However, I am still getting large errors in estimates of the gradient even
> when the mesh is nearly perfectly orthogonal. These errors DO NOT become
> smaller as resolution is increased. I have two question:
>
> Are these errors to be expected in unstructured meshes?

Cell centered FV has two issues with unstructured meshes. You know
about the orthogonality issue, but there is a second issue. The
docstring from the following code might help explain,


https://github.com/usnistgov/fipy/blob/develop/fipy/variables/cellVariable.py#L257

See, `\frac{1}{V_P} \sum_f \vec{n} \phi_f A_f`. This is the
discretized gradient calculation. In this formulation using the
triangular mesh, the calculated \phi_f is not the correct average for
the face as the line between the cell centers doesn't cross the face
in the middle of the face. I think this is called non-conjunctionality
in finite volume speak. See equation 9.8 in the following link for
further explanation (we should probably implement this approach).

   http://bit.ly/2aqRiHt

I believe that this is the source of the error in the gradient
calculation. It's a consistent error independent of the mesh
refinement I think (or maybe first order). An alternative gradient
calculation is the least squares gradient, see

   
https://github.com/usnistgov/fipy/blob/develop/fipy/variables/cellVariable.py#L270

Using that


import fipy as fp
nx = 25
dx = 1. / nx
mesh = fp.Tri2D(nx=nx, ny=nx, dx=dx, dy=dx)
psi = fp.CellVariable(mesh=mesh)
psi.constrain(1.0, where=mesh.facesLeft)
psi.faceGrad.constrain(1.0, where=mesh.facesRight)
eq = fp.DiffusionTerm().solve(psi)
print psi.leastSquaresGrad[:, :100]


seems to give lots of 1s and 0s which seems correct. The cells next to
the boundary are not correct as they aren't using the boundary
conditions in the calculation. I'm not sure about the order of
accuracy of the least squares or why it happens to work for this mesh.

> Is it acceptable to average the gradient over space to reduce the error?

I don't think so unless it can be shown to improve the accuracy with analysis.

-- 
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

hi,

i scanned the list and the documentation but have not been able to make
this work properly yet:

i have a 2D pde of the form

\[

\partial_t \rho =  \nabla \cdot (\rho \nabla \phi) + D \nabla^2 \rho +
\lambda \rho

\]

where phi is a potential field and lambda is a source field, both are
constant in time and smoothly varying in space.

i have been solving this with fipy, using

DiffusionTermCorrection(coeff=[D,]) for the diffusion
PowerLawConvectionTerm(coeff=v ) where v = -phi.grad for the advection
ImplicitSourceTerm(coeff=lambda) for the source

the boundary conditions are either no-flux or dirichlet. i think
advection is dominant, although i can't say i have actually tried to
quantify that.

i am testing the solution by subtracting the rhs and lhs of the equation
and comparing that to \partial_t \rho in magnitude. i am having trouble
making this error smaller than around 10%: at points where the density
grows strongly, the solution deviates.

at the moment i am using a loop like this:

for (st, stnext) in zip(Solve.sample_t, Solve.sample_t[1:]):
for tt in nx.arange(st, stnext, Solve.step):
rho.updateOld()
res = 1.
while res > Solve.tolerance:
res = eq.sweep(var=rho, dt=Solve.step)

i've played around with making the mesh smaller -- did not seem to make
any difference; reducing the tolerance -- some effect, and reducing the
temporal step size, which did not help much.

should i be doing something differently?

thanks for any help!

n.
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]