Material property depending on solution

1

Hi all,

I’d like to update material properties as function of a solution field (say temperature) on a several subdomains. So far, the only way I managed to do it is by iterating through the mesh and assigning values cell by cell.

from fenics import *


def coeff(T, subdomain_id):
    return [2*T, 3*T][subdomain_id]

nx = ny = 10
mesh = UnitSquareMesh(nx, ny)
n = FacetNormal(mesh)
V = FunctionSpace(mesh, 'P', 1)
D0 = FunctionSpace(mesh, 'DG', 0)

# Temperature
temp = Expression("2*x[0]*x[0] + t", degree=2, t=1)
# Solution
u = interpolate(temp, V)
# Material temperature dependent property
D = Function(D0)

# Time stepping loop
for t in [1, 2, 3]:
    temp.t = t

    # Update coefficient D
    for cell in cells(mesh):
        x = cell.midpoint()
        if x.x() > 0.5:
            subdomain_id = 1
        else:
            subdomain_id = 0
        cell_no = cell.index()
        D.vector()[cell_no] = coeff(temp(x), subdomain_id)

    # Compute flux
    print(assemble(D*dot(grad(u), n)*ds))

In this exemple, temp is the temperature on the domain, D is the property (say diffusion coefficient) that varies with temperature as shown in coeff().

This is really heavy especially for big meshes and I am quite sure there’s a clever way to do this.

Does anyone have an idea ?

Cheers
Rem

2

Hi Rem,

I’m really not an expert but have you tried a UserExpression?
Something like

class CoeffClass(UserExpression):
    def __init__(self, temp, **kwargs):
        super().__init__(kwargs)
        self.temp = temp
    def eval(self, values, x):
        values[0] = ... 2*self.temp(x) ...
        if x[0] > 0.5:
            values[1] = 1
        else:
            ...
            values[1] = 0
    def value_shape(self):
        return ()

Indexing needs to match your function space dimensions, of course.
Then

D = CoeffClass(temp=temp)
...
# Update coefficient D
D.temp = temp

The update step might not even be necessary.
Conditional statements from the UFL Package might even be a much better/faster solution but I never used them in combination with spatial coordinates and your kinds of function spaces.

Hope it helps at least a little :sweat_smile:

5

Thanks ! it works perfectly ! here’s the MWE

from fenics import *

class CoeffClass(UserExpression):

    def __init__(self, temp, **kwargs):

        super().__init__(kwargs)

        self.temp = temp

    def eval(self, value, x):

        if x[0] > 0.5:

            value[0] = 2*self.temp(x)

        else:

            value[0] = 3*self.temp(x)

    def value_shape(self):

        return ()

nx = ny = 100

mesh = UnitSquareMesh(nx, ny)

n = FacetNormal(mesh)

V = FunctionSpace(mesh, 'P', 1)

D0 = FunctionSpace(mesh, 'DG', 0)

# Temperature

temp = Expression("2*x[0]*x[0] + t", degree=2, t=1)

D = CoeffClass(temp=temp, degree=1)

# Solution

u = interpolate(temp, V)

# Time stepping loop

for t in [1, 2, 3]:

    temp.t = t

    # Update coefficient D

    D.temp = temp

    # Compute flux

    print(assemble(D*dot(grad(u), n)*ds))
1 Like