How to add a ceiling (limit value) on 1D Diffusion case.

[beantsja]

Hi, I am trying to solve 1D Diffusion equation, where there is a limit of the concentration value over the domain, as if there is a ceiling.

This picture is basically 50% of what I am trying to achieve. Because, in this picture, if the value at the index exceeds the limit value, the limit value simply replaces the value at the index using Constrain().

What I want to achieve is even when the value at the index exceeds the limit, these extra flux flow into the next domain, basically solving 1D diffusion but as if there is a ceiling that is limiting the concentration to grow (The top part is basically closed).

I was looking at the examples and not sure how to add this implementation.

[guyer]

What I'm used to seeing for a conserved quantity with limiting values is a diffusion coefficient that depends on the quantity, e.g.

$$ \frac{\partial C}{\partial t} = \nabla\cdot\left[ D_0 C (1 - C) \nabla C \right] $$

This keeps $C$ from going below 0 or above 1. In FiPy, this can be written as

eq = fp.TransientTerm(var=C) == fp.DiffusionTerm(coeff=D0 * C * (1 - C), var=C)

but be aware of the distinction between CellVariable and FaceVariable coefficients. You may want to write

eq = fp.TransientTerm(var=C) == fp.DiffusionTerm(coeff=D0 * C.faceValue * (1 - C).faceValue, var=C)

or even

eq = fp.TransientTerm(var=C) == fp.DiffusionTerm(coeff=D0 * C.harmonicFaceValue * (1 - C).harmonicFaceValue, var=C)

(C.faceValue is a synonym for C.arithmeticFaceValue. C.geometricFaceValue is not helpful in my experience, but it's available).

[beantsja]

eqI = TransientTerm() == DiffusionTerm(coeff=D_H)

This is what I used for 1D Diffusion equation, where I don't have "fp.". What is the "fp." for?

Also, I initialized the initial condition as

phi = CellVariable(name="phi", mesh=mesh, value=float(initial_phi)). Can I still initialize with CellVariable, not FaceVariable?

[guyer]

import fipy as fp

[guyer]

Your solution variable must be a CelVariable. Your diffusion coefficient should be a FaceVariable, but if you give a CellVariable, FiPy will interpolate it to a FaceVariable.