Unable to interpolate function into function space

Hi everyone,

I am currently trying to run some code in parallel via mpirun but when I do, I obtain this error message

Traceback (most recent call last):
  File "target3Dmsh.py", line 69, in <module>
    T0 = interpolate(initial_T, Space)
  File "/usr/lib/python3/dist-packages/dolfin/fem/interpolation.py", line 71, in interpolate
    Pv.interpolate(v._cpp_object)
  File "/usr/lib/python3/dist-packages/dolfin/function/function.py", line 365, in interpolate
    self._cpp_object.interpolate(u)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to interpolate function into function space.
*** Reason:  Wrong size of vector.
*** Where:   This error was encountered inside FunctionSpace.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[18617,1],1]
  Exit code:    1
--------------------------------------------------------------------------

When I run this code normally, without mpirun, my code runs but takes a very long time to iterate through to the next time step due to the scale of the problem. How can this error be fixed to allow for fast processing of my script?

Also, I do not think this is an issue with RAM or memory limitations as I am currently running this code within a large server. I am using Fenics 2019.1.0

My full code is shown below

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
#from mshr import *
from ufl import as_tensor
from ufl import Index
import math
#import mshr

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(20)

# SOLVE THE HEAT TRANSFER PROBLEM FOR A 3D TARGET

#----------------
t_start = 0.0
t_end = 20
nstep = 20
dtime = (t_end - t_start)/ nstep  #time step
#-----------------
#Configure for dimensions of system
#----------THERMAL--PROPERTIES--[Element: W]--------
kappa = 0.0172         #conductivity [W/mm K] 
c = 0.132         #Specific Heat Capacity [J/gK]
rho = 0.01925         #Density [g/mm^3]
const = kappa /(c * rho)
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag
alpha = 4.3943e+4 #absorption coefficient (mm^-1)
#---------------------------------------


pi = 3.141592653589793
T_am = 298 #ambient vacuum temperature (K)
T_a4 = T_am**4
epsilon = 0.25  # material emissivity
sigma = 5.67E-14 # W/(mm**2.K**4)
es = epsilon*sigma
Pmax = 500 #peak power of laser in W
w = 0.5     #mm
R = 1.5    #mm
area = pi*w*w  #area of target (mm^2)
depth = 10/alpha #thickness of target (mm)
volume = area * depth #volume of target
Iden = Pmax / (pi*w*w*volume*rho)    #Intensity per unit mass
absorb_depth = 1/alpha

Laser = Expression('2*A*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w*w))',degree=3, A=0.2, alpha=alpha, Iden=Iden,w= 0.56) #power (w2 localises the z-coordinates)
#A =! emissivity for system in thermal equillibrium


# Create mesh
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)    
    
Space = FunctionSpace(mesh, 'P', 1)      #define finite element function space, defined via basis functions
VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
cells = MeshFunction('size_t',mesh,mesh.topology().dim()) #codimension of 0
facets = MeshFunction('size_t',mesh,mesh.topology().dim()-1) #codimension of 1
da = Measure('ds', domain=mesh, subdomain_data = facets, metadata = {'quadrature_degree': 2})  #area element
dv = Measure('dx', domain=mesh, subdomain_data = cells, metadata = {'quadrature_degree': 2})   #volume element

initial_T = Expression("Tini", Tini=T_am, degree=3) # extrapolate an expression for the temperature before heating
T0 = interpolate(initial_T, Space)
T = Function(Space)         # Temperature
V = TestFunction(Space)     # Test Function used for FEA
dT = TrialFunction(Space)   # Temperature Derivative
q0 = Function(VectorSpace)  # heat flux at previous time step
i = Index()
G = as_tensor(T.dx(i), (i))  #gradient of T
G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step 

q = as_tensor(dtime/(dtime + tau_q) * (tau_q/dtime*q0[i] - kappa*(1+tau_T/dtime)*G[i] + kappa*tau_T/dtime*G0[i]),(i)) #heat
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i) - rho*Laser*V ) * dv + es*(T**4 - T_a4)*V*da   #final form to solve
Gain = derivative(F, T, dT)    # Gain will be usedas the Jacobian required to determine the evolution of a dynamic system 

file_T = File('target3D/solution.pvd')
for t in np.arange(t_start,t_end,dtime):
	print( "Time", t)
	solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "umfpack", "relative_tolerance": 1e-4} }, form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
	file_T << (T,t)
	q_tmp = project(q, VectorSpace, solver_type= "superlu")
	q0.assign(q_tmp)
	T0.assign(T)      #change so that T0 is equal to T for the next time step

This code also uses a mesh given below:

R = 1.5;
lc  = .2;
Point(1) = {0,0,0,lc};
Point(2) = {R,0,0,lc};
Point(3) = {0,R,0,lc};
Point(4) = {-R,0,0,lc};
Point(5) = {0,-R,0,lc};

Circle(1) = {2,1,3};
Circle(2) = {3,1,4};
Circle(3) = {4,1,5};
Circle(4) = {5,1,2};

Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};

Extrude {0,0,-8} {
Surface{6};
}

Field[1] = Distance;
Field[1].NodesList = {1};

Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = lc / 10;
Field[2].LcMax = lc;
Field[2].DistMin = 0.1;
Field[2].DistMax = 4;

Field[3] = Min;
Field[3].FieldsList = {2};
Background Field = 3;

Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature = 0;

Thanks in advance

I believe this has been fixed in the development version of fenics, available at quay.io/fenicsproject/dev:latest