I'm trying to use FiPy to solve for the transient behaviour of a simple 2D particle diffusion scenario with a small circular sink (particle concentration = n). I would like the sink to be saturable, i.e., there should be a maximum possible loss rate inside the sink. At the same time, the concentration should not be allowed to go negative. The precise behaviour of the sink should be:
- When the total rate of particles entering the sink <= maximum rate, the sink absorbs all the particles inside it. This is like a Dirichlet boundary condition at the sink with n = 0.
- When the total rate of particles entering the sink > maximum rate, the sink only absorbs particles at the maximum rate. I believe this could be expressed as a Neumann (flux) boundary condition at the perimeter of the sink.
The maximum rate is expressed in units of s^-1, i.e., it is the loss rate of particles inside the sink (not a rate of concentration loss). What is the most computationally-efficient way to express this in FiPy?
Background
In this scenario there is an initial Gaussian distribution in n, whose subsequent behaviour is governed by:
eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant)
I'm dealing with orthotropic diffusion, so the diff_coeff variable is given in the form [[[D1, 0], [0, D2]]].
The aim is to see the impact of one small circular sink on this solution. Since the sink is much smaller than the initial distribution and diffusion length, I'm using a Gmsh mesh to reduce calculation time. In this mesh, the cell size becomes much smaller near the sink.
I've already done this for an "absolute" sink (i.e., one which fixes n to 0 locally by using a Dirichlet boundary condition on the perimeter of the sink):
n.constrain(value=0, where=mesh.physicalFaces["sink"])
Solving for n(t) in this case works very well, and I can take large implicit time steps without sacrificing too much accuracy.
My issue arises when I try to make the sink "zero-order", i.e., include it in the equation as:
sink = CellVariable(mesh=mesh)
sink.setValue(value=maximum_rate, where=sink_coords <= sink_radius)
eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant) - sink
Solving for n(t) in this case also works well with large implicit steps, but when the maximum_rate is large, or the sink is far from the initial distribution, n can be negative around the sink. For my scenario, n is not allowed to be negative, and this is the main problem I'm facing.
If I try to control things so n does not go negative, I am unable to take large implicit time steps without losing a lot of accuracy (see Full code example at the end). Is there a way to implement a saturable sink while still being able to take large timesteps with reasonable accuracy? (Similar sized timesteps to the "absolute" sink case).
Attempted solutions
So far I have only tried naïve ways to solve the problem, all of which have resulted in the necessity of much smaller timesteps in order to get a reasonably accurate solution. One is to simply force n to be non-negative after each timestep:
eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant) - sink
for step in range(steps):
eq.solve(var=n, dt=time_step)
n.setValue(value=0, where=n < 0)
The other is to completely exclude the sink from the FiPy equation, and instead manually account for it in each timestep:
eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant)
for step in range(steps):
eq.solve(var=n, dt=time_step)
n.setValue(value = n - maximum_rate * time_step, where=sink_coords <= sink_radius)
n.setValue(value=0, where=n < 0)
It isn't surprising that these methods need much smaller time_step. The first method overestimates the impact of the sink if time_step is large since (as far as I understand) in the eq.solve step diffusion can still be caused by negative n. The second method underestimates the sink impact with large time_step.
I imagine there should be some solution in which the dependence of sink on n could be explicitly expressed. Something like:
expression = n / time_step if n / time_step < maximum_rate else maximum_rate
sink.setValue(value=expression, where=sink_coords <= sink_radius)
eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant) - sink
But I can't think of a way to express it properly. I suppose I would need some numerix function which is like a Heaviside function with a constant slope in the transition region?
Another solution could also be in using flux boundary conditions instead, but I am unsure about how to proceed with this. Any help provided is much appreciated, and thanks in advance! What follows is one full example of what I am currently doing.
Full code example
Using Python 3.10.4, FiPy 3.4.3, Gmsh 4.6.0.
My current simple test case is comparing an "absolute" sink (one which absorbs all particles, using Dirichlet condition n = 0) with a saturable sink with a very high maximum_rate. Ideally, these scenarios should give identical solutions, since the sink should behave as an absolute one until the total rate of particles entering the sink > maximum_rate - if maximum_rate is very high, this condition will never be met, so the sink should always be an absolute one.
from fipy import CellVariable, Gmsh2D, TransientTerm, DiffusionTerm, ImplicitSourceTerm
from fipy.tools import numerix
# Define some parameters.
Ds = [50000, 12500] # diffusion constants
sink_r = 1 # sink radius
initial_fwhm = 70 # initial state FWHM
max_time = 2 # up to which time we care about
x_sink = 300 # sink x coordinate
y_sink = 0 # sink y coordinate
constant = 4 # Background decay
use_dirichlet = False # if set to True, use n = 0 Dirichlet perimeter BC, no saturation.
# For testing purpose, set the maximum loss rate very high, so solution should match
# the use_dirichlet scenario.
maximum_rate = 1e5
# Define Gaussian function for initial condition.
def gauss(x, y, fwhm):
sigma = fwhm / 2 * numerix.sqrt(2 * numerix.log(2))
return numerix.exp(-(x**2 + y**2) / (2 * sigma**2)) / (sigma**2 * 2 * numerix.pi)
# Generate Gmsh.
major_radius = 2.5 * 2 * numerix.sqrt(Ds[0] * max_time) + 1.5 * initial_fwhm + 1.1 * sink_r
minor_radius = 2.5 * 2 * numerix.sqrt(Ds[1] * max_time) + 1.5 * initial_fwhm + 1.1 * sink_r
cell_size = minor_radius / 10
sink_cell_size = sink_r / 2
sink_extra_radius = sink_r * 3
sink_limit_radius = sink_extra_radius + 3 * (cell_size - sink_cell_size)
initial_cell_size = initial_fwhm / 4
initial_limit_radius = minor_radius * 0.66
# Set the sink as a "hole" in the surface if we are using n=0 BCs.
plane_surface = "19, 20" if use_dirichlet else "19"
embed = "" if use_dirichlet else "15, 16, 17, 18"
gmsh_script = f"""
Point(1) = {{0, 0, 0, {initial_cell_size}}};
Point(2) = {{0, {-minor_radius}, 0, {cell_size}}};
Point(3) = {{{major_radius}, 0, 0, {cell_size}}};
Point(4) = {{0, {minor_radius}, 0, {cell_size}}};
Point(5) = {{{-major_radius}, 0, 0, {cell_size}}};
Ellipse(6) = {{2, 1, 3}};
Ellipse(7) = {{3, 1, 4}};
Ellipse(8) = {{4, 1, 5}};
Ellipse(9) = {{5, 1, 2}};
Point(10) = {{{x_sink}, {y_sink}, 0, {sink_cell_size}}};
Point(11) = {{{-sink_r + x_sink}, {y_sink}, 0, {sink_cell_size}}};
Point(12) = {{{x_sink}, {sink_r + y_sink}, 0, {sink_cell_size}}};
Point(13) = {{{sink_r + x_sink}, {y_sink}, 0, {sink_cell_size}}};
Point(14) = {{{x_sink}, {-sink_r + y_sink}, 0, {sink_cell_size}}};
Circle(15) = {{11, 10, 12}};
Circle(16) = {{12, 10, 13}};
Circle(17) = {{13, 10, 14}};
Circle(18) = {{14, 10, 11}};
Line Loop(19) = {{6, 7, 8, 9}};
Line Loop(20) = {{15, 16, 17, 18}};
Plane Surface(21) = {{{plane_surface}}};
Curve{{{embed}}} In Surface{{21}};
Physical Surface("main") = {{21}};
Physical Line("sink") = {{15, 16, 17, 18}};
Field[1] = Distance;
Field[1].NodesList = {{1}};
Field[2] = Distance;
Field[2].NodesList = {{10}};
Field[3] = Threshold;
Field[3].IField = 1;
Field[3].LcMin = {initial_cell_size};
Field[3].LcMax = {cell_size};
Field[3].DistMin = 0;
Field[3].DistMax = {initial_limit_radius};
Field[4] = Threshold;
Field[4].IField = 2;
Field[4].LcMin = {sink_cell_size};
Field[4].LcMax = {cell_size};
Field[4].DistMin = {sink_extra_radius};
Field[4].DistMax = {sink_limit_radius};
Field[5] = Min;
Field[5].FieldsList = {{3, 4}};
Background Field = 5;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature = 0;
"""
mesh = Gmsh2D(gmsh_script)
# Set up n.
n = CellVariable(name="solution variable", mesh=mesh, value=0.0)
x, y = mesh.cellCenters
n.setValue(gauss(x, y, 70))
# Set up sink.
sink_coords = numerix.sqrt((x - x_sink) ** 2 + y**2)
adj_rate = maximum_rate / (numerix.pi * sink_r**2) # Adjust for the size of the sink
sink = CellVariable(mesh=mesh)
sink.setValue(value=adj_rate, where=sink_coords <= sink_r)
# Set up equation.
diff_coeff = [[[Ds[0], 0], [0, Ds[1]]]]
eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant) - sink
if use_dirichlet: n.constrain(value=0, where=mesh.physicalFaces["sink"])
# Define time step based on mesh.
time_step = 0.25 * numerix.mean(mesh.cellVolumes) / (4 * min(Ds))
steps = int(max_time / time_step)
# Solve up to max_time.
for step in range(steps):
eq.solve(var=n, dt=time_step)
n.setValue(value=0, where=n<0)
If the n returned by this code with use_dirichlet = False (attempted saturable sink case) is compared to the n returned when use_dirichlet = True (absolute sink with n=0 boundary condition), the discrepancy is clear. Again, ideally these two cases should return the same solution for n at time max_time since the maximum_rate is set to be very high.
