title: Representing "zero-order" saturable sinks in FiPy efficiently and without negative concentration
body: <p>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 = <em>n</em>). I would like <strong>the sink to be saturable</strong>, i.e., there should be a maximum possible loss rate inside the sink. At the same time, <strong>the concentration should not be allowed to go negative</strong>. The precise behaviour of the sink should be:</p>
  <ul>
  <li>When the total rate of particles entering the sink &lt;= maximum rate, the sink absorbs all the particles inside it. This is like a Dirichlet boundary condition at the sink with <em>n</em> = 0.</li>
  <li>When the total rate of particles entering the sink &gt; 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.</li>
  </ul>
  <p>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). <strong>What is the most computationally-efficient way to express this in FiPy?</strong></p>
  <h2>Background</h2>
  <p>In this scenario there is an initial Gaussian distribution in <em>n</em>, whose subsequent behaviour is governed by:</p>
  <p><code>eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant)</code></p>
  <p>I'm dealing with orthotropic diffusion, so the <code>diff_coeff</code> variable is given in the form <code>[[[D1, 0], [0, D2]]]</code>.</p>
  <p>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.</p>
  <p>I've already done this for an &quot;absolute&quot; sink (i.e., one which fixes <em>n</em> to 0 locally by using a Dirichlet boundary condition on the perimeter of the sink):</p>
  <p><code>n.constrain(value=0, where=mesh.physicalFaces[&quot;sink&quot;])</code></p>
  <p>Solving for <em>n</em>(<em>t</em>) in this case works very well, and I can take large implicit time steps without sacrificing too much accuracy.</p>
  <p>My issue arises when I try to make the sink &quot;zero-order&quot;, i.e., include it in the equation as:</p>
  <pre><code>sink = CellVariable(mesh=mesh)
  sink.setValue(value=maximum_rate, where=sink_coords &lt;= sink_radius)
  eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant) - sink
  </code></pre>
  <p>Solving for <em>n</em>(<em>t</em>) in this case also works well with large implicit steps, but when the <code>maximum_rate</code> is large, or the sink is far from the initial distribution, <em>n</em> can be negative around the sink. <strong>For my scenario, <em>n</em> is not allowed to be negative</strong>, and this is the main problem I'm facing.</p>
  <p>If I try to control things so <em>n</em> does not go negative, I am unable to take large implicit time steps without losing a lot of accuracy <strong>(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?</strong> (Similar sized timesteps to the &quot;absolute&quot; sink case).</p>
  <h2>Attempted solutions</h2>
  <p>So far I have only tried naïve ways to solve the problem, all of which have resulted in the necessity of <em>much</em> smaller timesteps in order to get a reasonably accurate solution. One is to simply force <em>n</em> to be non-negative after each timestep:</p>
  <pre><code>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 &lt; 0)
  </code></pre>
  <p>The other is to completely exclude the sink from the FiPy equation, and instead manually account for it in each timestep:</p>
  <pre><code>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 &lt;= sink_radius)
      n.setValue(value=0, where=n &lt; 0)
  </code></pre>
  <p>It isn't surprising that these methods need much smaller <code>time_step</code>. The first method overestimates the impact of the sink if <code>time_step</code> is large since (as far as I understand) in the <code>eq.solve</code> step diffusion can still be caused by negative <em>n</em>. The second method underestimates the sink impact with large <code>time_step</code>.</p>
  <p>I imagine there should be some solution in which the dependence of <code>sink</code> on <code>n</code> could be explicitly expressed. Something like:</p>
  <pre><code>expression = n / time_step if n / time_step &lt; maximum_rate else maximum_rate
  sink.setValue(value=expression, where=sink_coords &lt;= sink_radius)
  eq = TransientTerm() == DiffusionTerm(diff_coeff) - ImplicitSourceTerm(1/constant) - sink
  </code></pre>
  <p>But I can't think of a way to express it properly. I suppose I would need some <code>numerix</code> function which is like a Heaviside function with a constant slope in the transition region?</p>
  <p>Another solution could also be in using flux boundary conditions instead, but I am unsure about how to proceed with this. <strong>Any help provided is much appreciated, and thanks in advance!</strong> What follows is one full example of what I am currently doing.</p>
  <h2>Full code example</h2>
  <p>Using Python 3.10.4, FiPy 3.4.3, Gmsh 4.6.0.</p>
  <p>My current simple test case is comparing an &quot;absolute&quot; sink (one which absorbs all particles, using Dirichlet condition <em>n</em> = 0) with a saturable sink with a very high <code>maximum_rate</code>. 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 &gt; <code>maximum_rate</code> - if <code>maximum_rate</code> is very high, this condition will never be met, so the sink should always be an absolute one.</p>
  <pre><code>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 &quot;hole&quot; in the surface if we are using n=0 BCs.
  plane_surface = &quot;19, 20&quot; if use_dirichlet else &quot;19&quot;
  embed = &quot;&quot; if use_dirichlet else &quot;15, 16, 17, 18&quot;
  
  gmsh_script = f&quot;&quot;&quot;
      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(&quot;main&quot;) = {{21}};
      Physical Line(&quot;sink&quot;) = {{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;
      &quot;&quot;&quot;
  
  mesh = Gmsh2D(gmsh_script)
  
  # Set up n.
  n = CellVariable(name=&quot;solution variable&quot;, 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 &lt;= 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[&quot;sink&quot;])
  
  # 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&lt;0)
  </code></pre>
  <p>If the <code>n</code> returned by this code with <code>use_dirichlet = False</code> (attempted saturable sink case) is compared to the <code>n</code> returned when <code>use_dirichlet = True</code> (absolute sink with <em>n</em>=0 boundary condition), the discrepancy is clear. Again, ideally these two cases should return the same solution for <code>n</code> at time <code>max_time</code> since the <code>maximum_rate</code> is set to be very high.</p>