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:
- Thermo-elastic evolution problem (full coupling) — Numerical tours of continuum mechanics using FEniCS master documentation
and - 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?