Perlin Noice Permeability Expression

I/O
1

Hi all,

I have an issue and honestly I am not sure, if a solution exists…

I am trying to model a Darcy flow through porous media: -div( k/mu * grad p ). This is fine, I have a simple demo that works. Now I wish to create a more complicated expression for the porosity, i.e. create another formulation for the permeability tensor k.

In my working demo, k is defined as the expression:

k=max_value(Expression('exp(-pow((x[1]-0.5-0.1*sin(10*x[0]))/0.1, 2))',degree=2), Expression('0.01',degree=2))

The above formulation for k defines a porous fraction following a sine curve.

Now I wish to model a flow trough gravel. For this I am using Perlin Noise and the code from: https://stackoverflow.com/questions/42147776/producing-2d-perlin-noise-with-numpy?fbclid=IwAR0IUD9FyM9DbbwRXhg9KmdGOUvy275h1y-JYnXDApa-nWmQ5EUx29Oi91Y

The definition of the Perlin Noice can be used within FEniCS by help of an UserExpression:

class RandomField(UserExpression): 
    def eval(self, value, x):
        self.s=5
        value[:] = perlin(x[0]*self.s, x[1]*self.s, seed=8465)/2. + 0.5
    def value_shape(self):
        return ()

k=RandomField(degree=2)

My issue is, that every time k is used it is generated again, which is very time consuming. Can it be implemented such that k is only initialized once and then used as a tensor afterwards?

Thank you.

Below is a minimal working example, showing that the test-UserExpression is executed twice

 from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

#Create mesh
mesh = UnitSquareMesh(64, 64)
n = FacetNormal(mesh)
#create function space
Q = FunctionSpace(mesh, 'DG', 0) # pressure space, Discontinuous Lagrange space


class test(UserExpression):
    def eval(self, values, x):
        values[:] = np.random.rand()
    def value_shape(self):
        return ()

k=test(degree=2)
kplot = project(k,Q)
f1 = plt.figure(1)
c=plot(kplot)
plt.colorbar(c)

kplot = project(k,Q)
f1 = plt.figure(2)
c=plot(kplot)
plt.colorbar(c)

plt.show()
2

You want to define a permeability field that is a function of space. I suggest you project the expression to , for example, DG space

DG_space = FunctionSpace(mesh, 'DG', 0)
k = project(your_expression_object, DG_space)

Then you can use that k in your variational formulation without having to compile a different expression.

I hope this is helpful!

3

Thank you, that works perfectly!