Non-linear solver error for thermo-elastic equation

Hi Everyone

I’m coding a non-lineal thermo-mechanical model where the coefficient of thermal expansion and the Young modulus depend on the temperature. I used the demos:

  1. Thermo-elastic evolution problem (full coupling) — Numerical tours of continuum mechanics using FEniCS master documentation
    and
  2. https://fenicsproject.org/olddocs/dolfin/1.4.0/python/demo/documented/nonlinear-poisson/python/documentation.html

However I could not solve this error:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F’.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.849e+05 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Traceback (most recent call last):
File “non_linear_thermoelasticity_fenics.py”, line 144, in
solve(form == 0, U, bcs, solver_parameters={“newton_solver”:{“relative_tolerance”: 1e-6}})
File “/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/dolfin/fem/solving.py”, line 266, in _solve_varproblem
solver.solve()
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 successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics_1555852928859/work/dolfin/la/PETScMatrix.cpp.
*** Process: 0


*** DOLFIN version: 2018.1.0
*** Git changeset: 9e1777c0df7c762c7230189c87c13fd5ad33a4f0
*** -------------------------------------------------------------------------

Here is the full code and in this link you can download the mesh and other files to run it (Dropbox - Files.zip - Simplify your life) :

from dolfin import *
from mshr import *
from fenics import *
import numpy as np
if has_linear_algebra_backend("Epetra"):
    parameters["linear_algebra_backend"] = "Epetra"

has_linear_algebra_backend('PETSc')
# 1. Read mesh
mesh = Mesh("dyke_mesh_super_fine.xml")
facets = MeshFunction("size_t", mesh, "dyke_mesh_super_fine_facet_region.xml")
ds = ds(subdomain_data=facets)
#caracteristica de la malla
#W=6000 #Model width m
#H=3000 #Model height m
#Prof=1500 #Mean depth magma chamber m
#w1=10 #Horizontal length of the magma chamber m
#h1=500 #Vertical length of the magma chamber m
#el (0,0) del modelo está en el vertice superior izquierdo

#2. Read rock properties
data=np.loadtxt('Model_catalog.csv',skiprows=1,delimiter=',',usecols = (2,3,4,5,6,7,8) )
Temperatura=data[:,0]
Moduloelasticidad=data[:,1]
Poisson=data[:,2]
Expansiontermica=data[:,3]
Capacidadcalorifica=data[:,4]
Conductividadtermica=data[:,5]
Presionmagma=data[:,6]
nmodel=len(Temperatura)

for i in np.arange(nmodel):
	tmax=24*5   #horas
	num_step=10
	dtime=tmax/num_step
	Tintru=Temperatura[i]  #grados celsius
	rho0=2700 #densidad
	E0=Moduloelasticidad[i]  #modulo de elasticidad en MPa
	nu0=Poisson[i]  # razon de poisson, adimensional
	alpha0=Expansiontermica[i] #coeficiente de expansion termica, 1/°C
	cv0=Capacidadcalorifica[i] #conductividad termica, J/kg°C
	k0=Conductividadtermica[i] #conductividad térmica, J/horas/m°C
	p0=Presionmagma[i] #Magma pressure, MPa
	filename="sol_thermo_Tintro="+str(Tintru)+"_rho="+str(rho0)+"_E0="+str(E0)+"_nu0="+str(nu0)+"_alpha0="+str(alpha0)+"_cv0="+str(cv0)+"_k0="+str(k0)+"_p0="+str(p0)+".xdmf"
	Tm=Constant(Tintru) #temperatura magma
	T0=Expression('15-0.03*x[1]',degree=1) #gradiente termal corteza
	DTdike=Expression('Tm-T0',degree=1,Tm=Tm,T0=T0)
	E=Constant(E0)
	nu = Constant(nu0)
	lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
	mu = Constant(E/2/(1+nu))
	rho = Constant(rho0)     # density
	alpha = Constant(alpha0)  # 8e-6 thermal expansion coefficient
	#kappa = Constant(alpha*(2*mu + 3*lmbda))
	cV = Constant(cv0)*rho  # specific heat per unit volume at constant strain
	k = Constant(k0)  # thermal conductivity
	rho_g = 0   # 0 si no hay fuerzas de cuerpo, 0.027 si las hay
	f = Constant((0, -rho_g))
	#P0=Expression(("-p0*(x[0]-3000)/sqrt(pow(x[0]-3000,2)+pow(5/250,4)*pow(x[1]+1500,2))","-p0*pow(5/250,2)*(x[1]+1500)/sqrt(pow(x[0]-3000,2)+pow(5/250,4)*pow(x[1]+1500,2))"),degree=1,p0=p0)
	q=Constant(p0)

	def kappa(T):
		return alpha0*T*(2*mu + 3*lmbda)

	#4. Define space of solutions and boundary conditions
	Vue = VectorElement('CG', mesh.ufl_cell(), 2) # displacement finite element
	Vte = FiniteElement('CG', mesh.ufl_cell(), 1) # temperature finite element
	V = FunctionSpace(mesh, MixedElement([Vue, Vte]))

	bc1 = DirichletBC(V.sub(0), Constant((0.,0.)), facets,102) #inferior
	bc2 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets,101) #laterales
	bc4 = DirichletBC(V.sub(1), DTdike, facets,104) #interior
	bc5 = DirichletBC(V.sub(1), Constant(0.), facets,102) #inferior
	bc6 = DirichletBC(V.sub(1), Constant(0.), facets,103) #superficie
	bc7 = DirichletBC(V.sub(1), Constant(0.), facets,101) #laterales
	bcs = [bc1, bc2, bc4, bc5, bc6, bc7]
	n = FacetNormal(mesh)

	#5.Espacios de funciones y formulaciones debiles
	U_ = TestFunction(V)
	(u_, Theta_) = split(U_)
	dU = Function(V)
	#dU = TrialFunction(V)
	(du, dTheta) = split(dU)
	Uold = Function(V)
	(uold, Thetaold) = split(Uold)




	def eps(v):
	    return sym(grad(v))

	def sigma(v, Theta):
		return (lmbda*tr(eps(v)) - kappa(Theta)*Theta)*Identity(2) + 2*mu*eps(v)

	dt=Constant(dtime)
	T0=project(T0,V.sub(1).collapse())

	mech_form = inner(sigma(du, dTheta), eps(u_))*dx - inner(f, u_)*dx + q*dot(n,u_)*ds(104)

	therm_form = (cV*(dTheta-Thetaold)/dt*Theta_ +
		      kappa(dTheta)*T0*tr(eps(du-uold))/dt*Theta_ +
		      dot(k*grad(dTheta), grad(Theta_)))*dx
	form = mech_form + therm_form



	U = Function(V)
	file_results = XDMFFile(filename)
	file_results.parameters["flush_output"] = True
	file_results.parameters["functions_share_mesh"] = True
	t=0


	for i in np.arange(num_step):
		t=t+dtime
		solve(form == 0, U, bcs, solver_parameters={"newton_solver":{"relative_tolerance": 1e-6}})
		#J = derivative(form, U)
		#solve(form == 0, U, bcs, J=J)
		#R  = action(form,U)
		#DR  = derivative(form,U)   # Gateaux derivative
		#problem = NonlinearVariationalProblem(R, U, bcs, DR)
		#solver  = NonlinearVariationalSolver(problem)
		Uold.assign(U)
		#u, T = split(U)
		u, T = U.split()
		u.rename('Displacement','label')
		file_results.write(u, t)
		T.rename('Temperature variation','label')
		file_results.write(T,t)
		Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
		sig = Function(Vsig, name="Cauchy_Stress")
		sig.assign(project(-sigma(u,T), Vsig))
		file_results.write(sig, t)
		Epssig = TensorFunctionSpace(mesh, "DG",degree=1)
		epsilon = Function(Epssig, name="Strain")
		epsilon.assign(project(eps(u), Epssig))
		file_results.write(epsilon,t)
		sx = sig.sub(0)
		sy = sig.sub(3)
		txy = sig.sub(1)

		def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
				return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))

		def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
				return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))

		S1 = sigma1(sx,sy,txy)
		S3 = sigma3(sx,sy,txy)
		Difstress=S1-S3
		Promstress=(S1+S3)/2
		Vp = FunctionSpace(mesh, "P",1)
		sig1=Function(Vp,name="Sigma 1")
		sig1.assign(project(S1,Vp))
		sig3=Function(Vp,name="Sigma 3")
		sig3.assign(project(S3,Vp))
		sigdif=Function(Vp,name="Differential stress")
		sigdif.assign(project(Difstress,Vp))
		sigprom=Function(Vp,name="Mean Stress")
		sigprom.assign(project(Promstress,Vp))
		file_results.write(sig1,t)
		file_results.write(sig3,t)
		file_results.write(sigdif,t)
		file_results.write(sigprom,t)
		print('Progress: '+str(100*t/tmax)+'%')

Some idea to fix it?

Please make sure that the code is formatted with 3x` such as

```python
from dolfin import *
# add code here
```

Also, it is quite hard to give advice when the error cannot be reproduced (as you have not supplied the mesh, or used a built in mesh). Please check if you can reproduce the error on a UnitCubeMesh or a UnitSquareMesh

You have also not supplied the full error message.

Hi Dokken, now I think it is properly formatted. I also add the files (meshes and .csv necessary to run it)

Your form is not defined with U, but dU which should therefore be the input to solve