masked diffusion term, or gmsh physical groups

[Urzainqui Aramburu Iñaki (Luke)]

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

[Daniel Wheeler]

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

> --
> 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.



--
Daniel Wheeler

[Daniel Wheeler]

Sending to the fipy list...

On Wed, Dec 8, 2021 at 10:41 AM Daniel Wheeler
<daniel....@gmail.com> wrote:
>
> "repr(coeff.faceValue)" is showing you the dependency tree of the
> variable "coeff.faceValue". For example,
>
> ~~~~
> >>> from fipy import Variable
> >>> a = Variable([1, 2])
> >>> b = Variable([2, 3])
> >>> c = a * b
> >>> repr(c)
> '(Variable(value=array([1, 2])) * Variable(value=array([2, 3])))'
> >>> print(c)
> [2 6]
> >>> c.value
> array([2, 6])
> >>> c
> (Variable(value=array([1, 2])) * Variable(value=array([2, 3])))
> ~~~~
>
> So, if we ask for the representation of "c" then we see the variables
> that "c" depends on, not the evaluation of c. Printing "c" directly or
> asking for its value forces an evaluation.
>
> On Wed, Dec 8, 2021 at 8:21 AM Urzainqui Aramburu Iñaki (LUKE)
> <inaki.u...@luke.fi> wrote:
> >
> > Hi again,
> >
> > Thanks a lot for your support! This solved my issues.
> >
> > Small comment:
> > After a bit of tinkering, it seems like coeff.faceValue does indeed return a FaceVariable object. However, if I run print(repr(coeff.faceValue)), it does not look like it gives a lazy evaluation of the dependent variables. Instead, I get a bunch of values (numbers) printed out.
> >
> >
> > Iñaki
> > ________________________________
> > From: Daniel Wheeler <daniel....@gmail.com>
> > Sent: Tuesday, November 30, 2021 09:16 PM
> > To: Urzainqui Aramburu Iñaki (LUKE) <inaki.u...@luke.fi>
> > Cc: fi...@list.nist.gov <fi...@list.nist.gov>
> > Subject: Re: [fipy] masked diffusion term, or gmsh physical groups

> --
> Daniel Wheeler



--
Daniel Wheeler