fipy | FiPy is a Finite Volume PDE solver written in Python

This question is similar to cellvariable*Diffusion in fipy, but for convection instead of diffusion. The solution is the same though. Use,

to convert the red terms into a convection term with a coefficient of f and a source term.

For the blue terms do exactly the same to get a transient term with a coefficient and a source term.

Edit: if we assume that then we can still use

and approximate explicitly (with f as a variable in FiPy and use grad). However, we can take it further and use,

and approximate the final explicitly.

Again, we can go even further with,

Again the last term can be solved explicitly. Depending on the form of f the last them should become more insignificant and thus, the explicitness less of an issue.

That won't work. Diffusion is not a reversible operator. Try to define the problem with maths and then we can help to implement the equations in FiPy.

These:

Having both a Dirichlet and a Neumann on the same face in FV turns the nature of the problem from a boundary value problem into an initial value problem. In that sense the problem becomes over-specified as the right hand side boundary condition is still required. There may be ways to handle it with FD/FV with some hacks. However, FiPy certainly isn't set up to handle this type of problem.

In this line, change

Writing efficient JAX code is very similar to writing efficient NumPy code: generally if you are using a for loop over rows of your data, your code will not be very efficient. Instead, you should strive to write your computations in terms of vectorized operations.

In your code, it looks like you are relying on many non-JAX elements (e.g. NumPy masked arrays, operations in FiPy, etc.) so it's unlikely that JAX will be able to improve your runtime. I'd focus instead on rewriting your code to make efficient use of NumPy, replacing the for-loop logic with NumPy vectorized operations.

Here is an example of expressing your function in terms of vectorized operations:

FiPy doesn't provide any capability like this and isn't likely to. If your meshes were simply connected and convex, then you could use scipy.spatial.tsearch with the mesh vertices. Unfortunately, based on your previous questions, this approach won't work for you.

Many years ago, I did some work with semiconductor heterostructures, where I needed to account for different domain geometries built up during successive lithographic and electrodeposition steps. I used the Shapely package to define abstract geometrical domains and then converted those domains to a FiPy mesh with Gmsh. Shapely can be used to check point containment, even when the domains are concave or disconnected. The specifics were a bit different from what you're doing, but I believe it's still applicable.

This gist illustrates building multiple complex domains, checking point containment in them, and then generating a single FiPy mesh from them (which was my use-case; you could generate separate meshes from separate domains).

Here is a refinement on the brute-force solution I initially offered:

FiPy does not provide a mesh class for what you want. UniformGrid2D could be used as the starting point for an equilateral triangular mesh (or a hexahedral mesh, if really required). StackOverflow isn't the right venue for that, though. Please open an issue if you'd like to pursue that.

Gmsh will produce structured, triangular meshes, but I'm not sure it can be forced to produce regular, equilateral triangles.