Inhomogeneous Poisson equation with user defined expression

variational formulation
1

Hi guys, I am writing code for the inhomogeneous Poisson equation - \nabla \cdot K(x) \nabla u = f, where K(x) is a piecewise constant function (defined in the code below).

The issue I encountered is that the assembling of the weak form is extremely time-consuming. I believe it is due to the way I defined the expression for K(x). It’s probably like iterating over a NumPy array using a for loop, which is very slow in python.

Here is the minimal working code. I made a time comparison between running assembling (\nabla u, \nabla v) and (K(x) \nabla u, \nabla v ) in the code.

from fenics import *
import time
import numpy as np
import math

# Define Layer structure dielectric constant
dielectric_consts = [1,50,1,50,1,50]
class MyExpression(UserExpression):
    def eval(self, values, x): 
        num_layer = len(dielectric_consts)
        width = 1.0/num_layer

        if x[1] <= width: 
            values[0] = dielectric_consts[0]
        elif x[1] <= 2*width: 
            values[0] = dielectric_consts[1]
        elif x[1] <= 3*width: 
            values[0] = dielectric_consts[2]
        elif x[1] <= 4*width: 
            values[0] = dielectric_consts[3]
        elif x[1] <= 5*width: 
            values[0] = dielectric_consts[4]
        else: 
            values[0] = dielectric_consts[5]
    def value_shape(self):
        return() #scalar function  

mesh = UnitSquareMesh(400,400)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("100*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.01)", degree=2)

# Define weak form 1
K = MyExpression()
a_1 = inner(K*grad(u), grad(v))*dx # Assmbling is super slow with K, not optimized.
L_1 = f*v*dx

# Define weak form 2
a_2 = inner(grad(u), grad(v))*dx

begin = time.time()
A, b = assemble_system(a_1, L_1,None)
end = time.time() - begin
print("With K ... time taken: ",end)

begin = time.time()
A, b = assemble_system(a_2, L_1,None)
end = time.time() - begin
print("Without K ... time taken: ",end)

Any ideas on how to optimize this in FEniCS?

1 Like
2

This issue was addressed in the book Solving PDEs in Python
– The FEniCS Tutorial Volume 1, 4.3.3 Chaper 4.

The solution is to use C++ code and pass it to python expression. :slight_smile:

Update:
However, I am not able to compile to code after using C++ code in python. Still need help :joy: :joy: :joy:

image
image815×140 9.29 KB

To reproduce the error, please run the following code. Thank you!

from fenics import *

cppcode = """
class DielectricConstant : public Expression
{ 
  void eval(Array<double>& values, const Array<double>& x) const
  {
    int num_layer =  dielectric.size();
    double width = 1.0/num_layer;
    for (int i = 0; i < num_layer; i++){
        if( (x[1] >= i*width) && (x[1] <(i+1)*width) ){
          values[0] = dielectric[i];
          break;
        }
    }
  }
public: 
  std::vector<double> dielectric {1,50,1,50,1,50};
};
"""
mesh = UnitSquareMesh(200,200)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("100*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.01)", degree=1)

K_cpp = Expression(cpp_code = cppcode,degree = 1)

# Define weak form
a= inner(K_cpp*grad(u), grad(v))*dx

begin = time.time()
A, b = assemble_system(a, L_1,None)
end = time.time() - begin
print("With K cpp code ... time taken: ",end)
3

For the latest version of dolfin, I would suggest having a look at: Bitbucket