 # Perlin Noice Permeability Expression

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-0.5-0.1*sin(10*x))/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*self.s, x*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()``````

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.