Hi Ian,
One way to deal with the sharp discontinuity could be to have very
fine grid spacing near the discontinuity. The following code solves a
diffusion problem with a step in the diffusion coefficient. It only
uses 4 cells in 1D. 2 of the cells are very large and 2 are very
small. As the small cells get smaller, the accuracy increases,
~~~~
import fipy
import numpy
def run_problem(epsilon):
L = 1.
N = 4
dx = (L - epsilon * 2) / (N - 2)
dx_part = numpy.ones((N - 2) / 2) * dx
dx_array = numpy.concatenate([dx_part, [epsilon, epsilon], dx_part])
D1 = 1e-2
D2 = 1.
mesh = fipy.Grid1D(dx=dx_array)
var = fipy.CellVariable(mesh=mesh)
var.constrain(0, where=mesh.facesLeft)
var.constrain(1., where=mesh.facesRight)
diff_coeff = fipy.CellVariable(mesh=mesh)
diff_coeff[:] = D1
diff_coeff.setValue(D2, where=mesh.x > L / 2)
fipy.DiffusionTerm(diff_coeff).solve(var)
predicted_value = (var[N / 2 - 1] + var[N / 2]) / 2.
actual_value = D2 / (D1 + D2)
return abs(predicted_value - actual_value)
for epsilon in (1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9):
print epsilon, run_problem(epsilon)
~~~~
For me this gives:
0.1 0.00212603546249
0.01 0.000194041173913
0.001 1.92361668152e-05
0.0001 1.92195296744e-06
1e-05 1.92178367886e-07
1e-06 1.91937703509e-08
1e-07 1.54912405126e-09
1e-08 3.16000559053e-09
1e-09 2.76879312811e-08
There is a point of diminishing returns when cell volumes vary by many
orders of magnitude.
Cheers,
Daniel
On Fri, Sep 16, 2016 at 1:44 PM, Campbell, Ian
<[email protected]> wrote:
> Hello All,
>
>
>
> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard
> elliptic diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
>
>
>
> The diffusion coefficient, D, changes sharply at x = 1.0. Although this is
> quite similar to examples.diffusion.mesh1D, there is a key difference:
>
>
>
> We are seeking an instantaneous change in the coefficient value at that
> face, i.e. something of the effect that on one side of the face (i.e. from x
> = 0 up to and including x = 1.0) diffusion coefficient D has a value D1.
> Then there is an abrupt change in diffusion coefficient from x >= 1.0 (very
> notable, including the face at x = 1.0) and for the remainder of the mesh,
> diffusion coefficient D has a value D2. Currently, we have to either decide
> on whether the single face at x = 1.0 gets the value D1 or D2. However, when
> we run with either of these cases, the results are not accurate (when
> compared against a benchmark from a commercial FV solver). We are getting a
> smooth interpolation across the face , rather than the sharp change
> correctly predicted by the benchmarking code.
>
>
>
> With the code below, we have only been able to set the coefficient at that
> middle face to one of either D1 or D2. This doesn’t result in the desired
> instantaneous change in the coefficient value, but instead, the new
> coefficient value only applies from the NEXT face (which is a distance dx
> away) from this critical interior boundary.
>
>
>
> Coefficient = FaceVariable(mesh = 1D_mesh)
> # Define facevariable on mesh
>
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))
> # Set coefficient values in 1st half of mesh
>
> Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] <= 2.0))) # Set coefficient
> values in 2nd half of mesh
>
>
>
> An alternative with the inequalities adjusted, but that still only permits a
> single coefficient at that face:
>
> Coefficient = FaceVariable(mesh = 1D_mesh)
> # Define facevariable on mesh
>
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] < 1.0))
> # Set coefficient values in 1st half of mesh
>
> Coefficient.setValue(D2, where=((1.0 <= 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] <= 2.0))) # Set coefficient values
> in 2nd half of mesh
>
>
>
> Sincerely,
>
>
>
> - Ian
>
>
>
> P.S. Daniel, thank you very much for your previous reply on sweep speeds! It
> was very helpful.
>
>
>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
