Hi Iñaki,
Thanks for your questions. This shouldn't be a problem for FiPy.
There are some notable errors in the code you have in the image in
your email. Rather than point out the issues I'm just going to provide
an example that works.
~~~~
from fipy import CellVariable, Grid2D, FaceVariable, DiffusionTerm,
Viewer, TransientTerm
mesh = Grid2D(nx=2, ny=2, Lx=1.0, Ly=1.0)
phi = CellVariable(mesh=mesh, value=0.01, hasOld=True)
X, Y = mesh.faceCenters
coeff = phi.faceValue * (X != 0.5)
print(coeff)
# bcs on left side
phi.constrain(0.01, where=(X < 0.5) & (Y == 0.0))
phi.constrain(1, where=(X < 0.5) & (Y == 1.0))
# bcs on right side
phi.constrain(1, where=(X > 0.5) & (Y == 0.0))
phi.constrain(0.01, where=(X > 0.5) & (Y == 1.0))
eqn = (TransientTerm() == DiffusionTerm(coeff))
for i in range(10):
phi.updateOld()
eqn.solve(phi, dt=0.1)
Viewer(phi).plot()
print(phi)
input('stop')
~~~~
The internal faces in the Y direction have 0 diffusion in the above
while the diffusion coefficient is also dependent on phi. The left and
right side are therefore independent problems. Here the mask is X !=
0.5. In your example define the diffusion coefficient with
coeff = diff.faceValue * mask
diff.faceValue is a face variable. coeff will then depend on diff and
diff with depend on phi. You can check that your coeff has the correct
dependencies using "print(repr(coeff))" in your code to make
absolutely sure you trust the lazy evaluation. That will print the
dependency tree. Probably avoid having a function to construct the
diffusion coefficient. It shouldn't matter if done correctly, but
doesn't buy you anything as it won't be called more than once.
I hope that helps.
Cheers,
Daniel
On Thu, Nov 4, 2021 at 5:23 AM Urzainqui Aramburu Iñaki (Luke)
<inaki.urzain...@luke.fi> wrote:
>
> Dear FiPy developers,
>
> I am trying to solve a diffusion PDE with a non-linear diffusion term in a
> domain that contains 2 subdomains. Those subdomains are not specified in the
> mesh, i.e., the mesh treats everything as a single domain. I want to achieve
> diffusivity = 0 between the faces of one of the subdomains. In other words, I
> am trying to mask the diffusivity FaceVariable so that some faces (given, for
> instance, by the values of a fp.FaceVariable) are identically 0.
>
> My approach so far is the following:
>
> eq = fp.TransientTerm() == fp.DiffusionTerm(coeff=compute_diffusivity(phi))
>
> And compute_diffusivity is given by:
>
> def compute_diffusivity(phi):
> diff = fp.numerix.sqrt(2.0*phi) / fp.numerix.exp(4.0*phi**2)
> diff_face = fp.FaceVariable(mesh=mesh,
> value=diff.arithmeticFaceValue.value)
>
> # Create mask of subdomain, all ones except zeros where needed
> mask_face = fp.FaceVariable(...)
>
> # Apply mask to diffusivity and return.
> # Multiplication works because the target value at the mask is zero.
> return fp.FaceVariable(mesh=mesh, value=diff_face.value * mask_face.value)
>
> FiPy does solve the problem, but I am worried that it does not understand
> that the diffusion variable is non-linear. I suspect that the problem is that
> the .value method used to construct the FaceVariable is determining its value
> directly, without the lazy evaluation of phi.
>
> How can I solve this problem? What is the best way to set such a mask in FiPy?
>
> Another approach might be to define the subdomains within the mesh (using
> Gmsh 'physical groups', perhaps?) and then reading and using it in FiPy.
> However, I don't know how to do either of those things. (Extra info: I need
> my code to be useful both with meshes created within FiPy and by Gmsh).
>
> I would really appreciate your help. Thank you in advance!
>
> Iñaki
>
> --
> To unsubscribe from this group, send email to fipy+unsubscr...@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+unsubscr...@list.nist.gov.
--
Daniel Wheeler
--
To unsubscribe from this group, send email to fipy+unsubscr...@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+unsubscr...@list.nist.gov.
