Setting internal "constant flux" boundary conditions correctly

[Isaac Paten]

Hi, 

I am trying to solve a diffusion equation using FiPy, with a constant flux condition on internal boundaries - I have attached both a pdf describing the equation, and also a file with a MWE of my code. 

I am using a 2D geometry with two phases, labelled 0 and 1. Phase 1 is a square immersed in phase 0. I only need to solve the equation in phase 1, and the constant-flux BC exists only on the faces between cells labelled 0 and 1. 

I do not get any errors when I run the code, but the solution is 0 everywhere, so something is going wrong. The definition of my internal boundary seems ok (I plotted it). Is there something wrong with the way I am inputting the implicit source term? 

Thanks, 

Isaac 

[Guyer, Jonathan E. Dr. (Fed)]

Isaac -

What you've set up is an
[internal fixed gradient](https://pages.nist.gov/fipy/en/latest/USAGE.html#internal-fixed-gradient)
whereas your boundary condition defines a fixed flux.
I see that we don't actually cover this case, but it's actually easier to set up than fixed gradient IMO,
because the finite volume method naturally works with fluxes via the divergence theorem.

First, you'll want to turn off diffusion at your boundary:

```diff
D0 = 1
D = FaceVariable(mesh=mesh, value=D0)
+D.setValue(0., where=internal_boundary_mask)
```

Next , you need to specify the orientation of your flux.
`mesh.faceNormals` defines the orientation of faces *when pointing from cell 0 to cell 1*.
This is useful for FiPy internally, but cell 0 vs cell 1 doesn't mean much to humans.

Instead, we can take advantage of the knowledge you have about what's inside and outside
your domains of interest. By taking the gradient of your mask (after converting it from Boolean)
and normalizing it by the distance between adjacent cells, we obtain a unit vector that points in
(or out) of the domain of interest:

```diff
+solid_mask_var = CellVariable(mesh=mesh, value=solid_mask_flat * 1.)
+boundary_normals = solid_mask_var.faceGrad * mesh._cellDistances

+flux = 1.
+boundary_flux = internal_boundary_mask * flux * boundary_normals
```

Finally, add the divergence of this boundary flux as a source term:

```diff
-gradient = - 1 / D0
-
-largeValue = 1e10


eq = (TransientTerm() == DiffusionTerm(coeff=D)
- + DiffusionTerm(coeff=largeValue * internal_boundary_mask)
- - ImplicitSourceTerm((internal_boundary_mask * largeValue * gradient
- * mesh.faceNormals).divergence))
-
+ + boundary_flux.divergence)
```

For anything further, please start a GitHub Discussion: https://github.com/usnistgov/fipy/discussions
It's a much easier environment for talking about math and code than this mailing list.

________________________________________
From: fi...@list.nist.gov <fi...@list.nist.gov> on behalf of Isaac Paten <isaac...@gmail.com>
Sent: Saturday, May 3, 2025 07:49
To: FIPY
Subject: [fipy] Setting internal "constant flux" boundary conditions correctly

Hi,

Thanks,

Isaac

--
To unsubscribe from this group, send email to fipy+uns...@list.nist.gov

View this message at https://list.nist.gov/fipy
To unsubscribe from this group and stop receiving emails from it, send an email to fipy+uns...@list.nist.gov<mailto:fipy+uns...@list.nist.gov>.

[Isaac Paten]

Thanks very much! That works great for me.

Best, 

Isaac