Error in Mesh Refinement FEniCS 2019.1.0

General
1

Hi,

I’m trying to solve a coupled reaction-diffusion system where the reaction term depends on the position of a moving particle. I was attempting to modify my (working) program by refining the mesh at a given point each time-step.

Attempting to do so leads to the following error message -

------------------- Start compiler output ------------------------
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp: In member function 'virtual void ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578_cell_integral_main_otherwise::tabulate_tensor(double*, const double* const*, const double*, int) const':
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:98:18: error: redeclaration of 'const double J_c0'
   98 |     const double J_c0 = coordinate_dofs[0] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[2] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:91:18: note: 'const double J_c0' previously declared here
   91 |     const double J_c0 = coordinate_dofs[0] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[2] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:99:18: error: redeclaration of 'const double J_c3'
   99 |     const double J_c3 = coordinate_dofs[1] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[5] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:92:18: note: 'const double J_c3' previously declared here
   92 |     const double J_c3 = coordinate_dofs[1] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[5] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:100:18: error: redeclaration of 'const double J_c1'
  100 |     const double J_c1 = coordinate_dofs[0] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[4] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:93:18: note: 'const double J_c1' previously declared here
   93 |     const double J_c1 = coordinate_dofs[0] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[4] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:101:18: error: redeclaration of 'const double J_c2'
  101 |     const double J_c2 = coordinate_dofs[1] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[3] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~
/tmp/tmpehfx73wc/ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578.cpp:94:18: note: 'const double J_c2' previously declared here
   94 |     const double J_c2 = coordinate_dofs[1] * FE10_C0_D01_Q1[0][0][0] + coordinate_dofs[3] * FE10_C0_D01_Q1[0][0][1];
      |                  ^~~~

-------------------  End compiler output  ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/path/jitfailure-ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578
Traceback (most recent call last):
  File "/home/path/program.py", line 85, in <module>
    gradc = project(grad(c),V_vec)
            ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/projection.py", line 132, in project
    A, b = assemble_system(a, L, bcs=bcs,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/assembling.py", line 361, in assemble_system
    b_dolfin_form = _create_dolfin_form(b_form, form_compiler_parameters)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/assembling.py", line 56, in _create_dolfin_form
    return Form(form,
           ^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dijitso/jit.py", line 216, in jit
    raise DijitsoError("Dijitso JIT compilation failed, see '%s' for details"
dijitso.jit.DijitsoError: Dijitso JIT compilation failed, see '/home/path/jitfailure-ffc_form_d7859e3c92af28c07015f68bd9fe02b9dd242578' for details

The program I’m (attempting to run) is the following,

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

k = 10
t_ref = 1/k
D_c = 1e-9
l_ref = np.sqrt(D_c/k)
kappa=100 # 100
gamma_R = 1 # 0.1
Lx = 1
Ly = 1
T = 10*t_ref
Lambda = 4000
h = Lambda*gamma_R*l_ref**3/(kappa*t_ref**2)

dt = t_ref/100 # Time-step
Nt = int(round(T/float(dt)))

# Create mesh and define function space
nx = ny = 50
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition

particlePosition = np.zeros((Nt,2))
particlePosition[0,0] = 0.5
particlePosition[0,1] = 0.5

phi = np.zeros(Nt+1)
phi[0] = np.pi/3
speed = 0.005

pos_init = particlePosition[0,:]

sigma = 1e-2
c_0 = Expression('1/(sigma*pow(2*pi,1/2))*exp( -1/(2*pow(sigma,2))*(pow(x[0]-x0,2) + pow(x[1]-y0,2)) )',
degree=2,sigma=sigma,x0=pos_init[0],y0=pos_init[1])

c_D = Expression('0*x[0] + 0*x[1]',degree=1)
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, c_D, boundary)
# Define initial value
c_n = interpolate(c_0, V)
#u_n = project(u_D, V)
# Define variational problem
c = TrialFunction(V)
v = TestFunction(V)

a = (1 + k*dt)*c*v*dx + D_c*dt*dot( grad(c),grad(v) )*dx
L = c_n*v*dx

# Time-stepping
c = Function(V)
t = 0
grad_c = np.zeros((Nt,2))
mod_grad_c = np.zeros(Nt)
for n in range(Nt-1):
    pos_curr = particlePosition[n,:]
    pt = Point(pos_curr[0],pos_curr[1])
    
    markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)
    for cell in cells(mesh):
        markers[cell] = cell.contains(pt)

    mesh = refine(mesh,markers)

    # Compute solution
    [A,b] = assemble_system(a,L,bc)
    delta = PointSource(V,pt,1)
    delta.apply(b)
    solve(A, c.vector(), b)

    c_n.assign(c)

    V_vec = VectorFunctionSpace(mesh,"CG",1)
    gradc = project(grad(c),V_vec)
    grad_c[n,:] = gradc(( pos_curr[0],pos_curr[1] ))
    mod_grad_c[n] = np.sqrt(grad_c[n,0]**2 + grad_c[n,1]**2)
    dphi_dt = kappa/(gamma_R)*(-grad_c[n,0]*np.sin(phi[n])
                 + grad_c[n,1]*np.cos(phi[n]))
    phi[n+1] = dt*dphi_dt + phi[n]
        

    particlePosition[n+1,0] = speed*np.cos(phi[n+1])*dt + particlePosition[n,0]
    
    particlePosition[n+1,1] = speed*np.sin(phi[n+1])*dt + particlePosition[n,1]

plt.plot(particlePosition[:,0],particlePosition[:,1])

If the following lines are removed, the program runs fine

    markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)
    for cell in cells(mesh):
        markers[cell] = cell.contains(pt)

    mesh = refine(mesh,markers)

Any idea what’s causing the error? I used the same mesh-refinement routine on a Poisson-solver program and it worked fine.

2

I cannot reproduce your error. However, your mesh refinement is quite strange, as

doesn’t affect the solution of c, as the new mesh is a new variable not related with a and L. So it is unclear why you refine prior to computing the solution.

1 Like
3

I see - looks like I need to shore up my knowledge of FEniCS fundamentals. Thanks