Integrating over surface fails with JIT error (AssertionError expr.ufl_is_terminal_modifier)

variational formulation
1

I am using fenics version 2019.2.0.dev0 in python to solve a drift-diffusion PDE for current flow in semiconductors and have a solution for the PDE. I am trying to extract the current flow across a few internal and external boundaries by integrating the current density function over the boundary surface. I followed a few examples from the discourse to define such surface integrals however when trying to integrate the current density function the code fails with a
`` assert expr.ufl_is_terminal_modifier
AssertionError"

when performing an assemble. My code looks as follows

from dolfin import *
import numpy as np

# Define the mesh and function space
mesh = RectangleMesh(Point(0,0),Point(40,40),40,10) 
element = VectorElement('P',triangle,3,dim=3)
V = FunctionSpace(mesh, element)

# Define subdomains of interest
SubDomain1 = CompiledSubDomain( '(x[0]>= -1e-14)&&(x[0]<=20+1e-14)' )
SubDomain2 = CompiledSubDomain( '(x[0]>= 20-1e-14)&&(x[0]<=40+1e-14)' )

# Define the boundaries of interest
Boundary1 = CompiledSubDomain( 'near(x[0],0)' )
Boundary2 = CompiledSubDomain( 'near(x[0],20)' )
Boundary3 = CompiledSubDomain( 'near(x[0],40)' )

# Assign Id to boundaries
BoundaryId = MeshFunction('size_t',mesh,mesh.topology().dim()-1)

Boundary1.mark(BoundaryId,1)
Boundary2.mark(BoundaryId,2)
Boundary3.mark(BoundaryId,3)

# create measures to integrate over the marked boundaries (internal and external boundaries use ds and dS respectively)
dMeas_int = ds(domain = mesh,subdomain_data=BoundaryId)
dMeas_ext = dS(domain = mesh,subdomain_data=BoundaryId)

# Define the PDE to solve

############# Initial condition #####################
ic = Expression(('p0','n0','2.0-0.05*x[0]'),degree=3,p0=1e-18,n0=1e-4)
u0 = project(ic,V)
pi_p0, pi_n0, phi_0 = split(u0)

############# Dirichlet boundary conditions #########
DC1 = DirichletBC( V.sub(2),Constant(2.0),Boundary1 )
DC2 = DirichletBC( V.sub(2),Constant(0.0),Boundary3 )
############## Trial and solution functions ############
u = Function(V)
pi_p, pi_n, phi_i = split(u)
vi_p,vi_n, vi_phi = TestFunctions(V)
Rip = 0.
Rin=0.
gamma = 0.02585
dt = 0.01
############### Potential PDE ###################
ai_phi = (0.6454775 * dot(grad(phi_i),grad(vi_phi)) - (pi_p-pi_n)*vi_phi)*dx
Li_phi = 1e-4 * vi_phi * dx

######## Drift Diffusion Current ###############
Jip = gamma * grad(pi_p)+pi_p * grad(phi_i)
Jin = gamma * grad(pi_n)-pi_n * grad(phi_i)

######## Define the weak pde form #############
ai_p = (pi_p * vi_p + dot(Jip, grad(vi_p) ) * dt ) * dx
ai_n = (pi_n * vi_n + dot(Jin, grad(vi_n) ) * dt ) * dx
Lip = (pi_p0 * vi_p + Rip * vi_p * dt) * dx
Lin = (pi_n0 * vi_n + Rin * vi_n*dt) * dx
F = ai_p+ai_n+ai_phi - Lip - Lin- Li_phi

# Solve the PDE
solve(F==0,u,[DC1,DC2])

The code upto this point works fine and a solution to the problem is obtained in u.

# Integrate the current density to get total current flowing (the problematic part that failed)
    nx = Constant(np.array((1.,0.)))
    J = Jin-Jip
    Jx = dot(J('+'),nx)

    ########## Surface integral forms to get the total current #############
    I1 = Jx * dMeas_int(1)+Jx * dMeas_ext(1)    # integration over Boundary1  
    I2 = Jx * dMeas_int(2)+Jx * dMeas_ext(2)    # integration over Boundary2
    I3 = Jx * dMeas_int(3)+Jx * dMeas_ext(3)    # integration over Boundary3

    Curr1 = assemble(I1)
    Curr2 = assemble(I2)
    Curr3 = assemble(I3)

Here on the line of ``assemble(I1)" I get a very long run-time error saying

    Calling FFC just-in-time (JIT) compiler, this may take some time.
    Traceback (most recent call last):
    ...
    ...
    ...
    in balance_modified_terminal
        assert expr._ufl_is_terminal_modifier_
    AssertionError

The error gives no clue as to what is wrong with the code and I am stuck trying to evaluate the total current integral

Any idea what might be going wrong?

2

Hi, In the code above i dont see a definition of dMeas0 and dMeas1, only dMeas_int and dMeas_ext. Should the be the same?

3

Yes sorry, my bad. They should be the same. I updated the previous post to correct the discrepancy.

4

Your script is still not complete, as you are missing the following variables: gamma, dt, Rip, Rin, so your issue is not reproducible.

5

Sorry, added the corresponding terms and checked the minimal code to reproduce the error. It works now and reproduces the error.

6

Hi, another note is that you mixed up the internal and exterior measure. ds is exterior, dS is interior.
The problem with your code is that you are trying to restrict the integral over the outer boundary to one side. The problem is fixed by:

J = Jin-Jip
Jx = dot(J,nx)

I1 = Jx("+") * dMeas_int(1)+Jx * dMeas_ext(1)    # integrationver Boundary1  
I2 = Jx("+") * dMeas_int(2)+Jx * dMeas_ext(2)    # integration over Boundary2
I3 = Jx("+") * dMeas_int(3)+Jx * dMeas_ext(3)    # integration over Boundary3

However, note that the positive restriction is arbitrary, but consistent. See Integrating over an interior surface for a fix for this.

1 Like
7

Thanks, @dokken that solved the problem.