PETsc error code is : 63

Linear Algebra
1

I am trying to solve three pde’s with mixedElement. I suspect i did somethin wrong in applying boundary condition but i am new to Fenics so would really appreciate any kind of help.

from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
#%matplotlib notebook

L = 1.
R = 0.1
N = 50  # mesh density

domain = Rectangle(Point(0., 0.), Point(L, L)) - Circle(Point(0., 0.), R, 100)
mesh = generate_mesh(domain, N)

T0 = Constant(293.)
DThole = Constant(10.)
E = 70e3
nu = 0.3
lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu))
rho = Constant(2700.)     # density
alpha = 2.31e-5  # thermal expansion coefficient
kappa = Constant(alpha*(2*mu + 3*lmbda))
cV = Constant(910e-6)*rho  # specific heat per unit volume at constant strain
k = Constant(237e-6)  # thermal conductivity
T = 2.0            # final time
num_steps = 50     # number of time steps
dt = T / num_steps # time step size
A1 = 40
A2 = 0.5
E1 = 30000
E2 = 60000
m = 0.7
n = 1.2
Rg = 8.31
T0 = Constant(283.)

Vue = VectorElement('CG', mesh.ufl_cell(), 2) # displacement finite element
Vte = FiniteElement('CG', mesh.ufl_cell(), 1) # temperature finite element
Vce = FiniteElement('CG', mesh.ufl_cell(), 1) # Curing finite element
VT = FunctionSpace(mesh, MixedElement([Vue, Vte, Vce]))

def inner_boundary(x, on_boundary):
    return near(x[0]**2+x[1]**2, R**2, 1e-3) and on_boundary
def bottom(x, on_boundary):
    return near(x[1], 0) and on_boundary
def left(x, on_boundary):
    return near(x[0], 0) and on_boundary


bc1 = DirichletBC(VT.sub(0).sub(1), Constant(0.), inner_boundary)
bc2 = DirichletBC(VT.sub(0).sub(1), Constant(0.), bottom)
bc3 = DirichletBC(VT.sub(0).sub(1), Constant(0.), left)

bc4 = DirichletBC(VT.sub(1), Constant(1000.), inner_boundary)
bc5 = DirichletBC(VT.sub(1), Constant(400.), bottom)
bc6 = DirichletBC(VT.sub(1), Constant(400.), left)
bc7 = DirichletBC(VT.sub(2), Constant(1.), inner_boundary)
bc8 = DirichletBC(VT.sub(2), Constant(0.1), bottom)
bc9 = DirichletBC(VT.sub(2), Constant(0.1), left)
#bc4 = DirichletBC(VT.sub(1), Constant(100.), boundary)
bcT = [bc1, bc2, bc3, bc4, bc5, bc6, bc7, bc8, bc9]


u = Function(VT)
du, dTheta, dalpha = split(u)
v = TestFunction(VT)
u_, Theta_, Alpha_ = split(v)
Uold = Function(VT)
uold, Theta_old, Alpha_old = split(Uold)
f = Constant(5)
Theta = dTheta - T0
def eps(v):
    return sym(grad(v))


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


def K_1(dTheta):
    return A1*exp(-E1/(Rg*dTheta))
def K_2(dTheta):
    return A2*exp(-E2/(Rg*dTheta))
def S_c(dalpha, dTheta):
    return (K_1(dTheta) + (K_2(dTheta)*pow(dalpha, m)))*pow((1-dalpha),n)

dt = Constant(0.)

mech_form = inner(sigma(du, dTheta), eps(u_))*dx
therm_form = dTheta*Theta_*dx + dt*dot(k*grad(dTheta), grad(Theta_))*dx - Theta_old*Theta_*dx - S_c(Alpha_old, Theta_old)*Theta_*dt*dx
#therm_form = (cV*(dTheta-Theta_old)/dt*Theta_ + dot(k*grad(dTheta), grad(Theta_)))*dx #- Constant(50)*Theta_*dx
#cure_form = ((dalpha - Alpha_old)/dt)*Alpha_*dx - S_c(Alpha_old)*Alpha_*dx
cure_form = dalpha*Alpha_*dx - Alpha_old*Alpha_*dx - dt*S_c(Alpha_old, Theta_old)*Alpha_*dx
form = mech_form + therm_form + cure_form




#Pushp, Pushp_to_VT


Theta_o = Expression('293.', degree=1)
Alpha_o = Expression('0.0293', degree=2)
#uold, Theta_old, Alpha_old = Uold.split()




Theta_old_tmp = VT.sub(1).collapse()
Alpha_old_tmp = VT.sub(2).collapse()

Theta_old = interpolate(Theta_o, Theta_old_tmp)
Alpha_old = interpolate(Alpha_o, Alpha_old_tmp)
parameters['allow_extrapolation'] = True # TODO
Nincr = 100
t = np.logspace(1, 4, Nincr+1)
#Nx = 100
#x = np.linspace(-5, 10, Nx)
#T_res = np.zeros((Nx, Nincr+1))
C = Function(VT)
#u = Function(VT)
for (i, dti) in enumerate(np.diff(t)):
    print("Increment " + str(i+1))
    dt.assign(dti)
    solve(form == 0, C, bcT)
    print ("vdfv")
    
    Uold.assign(u)
    #T_res[:, i+1] = [C(xi, 0.)[1] for xi in x]
    #vtkfile << (T_res, t)

I am getting the following error:

Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
 *** Reason:  PETSc error code is: 63 (Argument out of range)"
2

Please edit the post so that the code is formatted properly, using “```” instead of the blockquote. Furthermore, the code must be complete, i.e. include all required imports.

3

Remove the definition of C here. The input to

should be the function you used to define your residual, i.e.

4

It worked. Thanks. but i was trying to follow the logic and procedure given in the following link which solves coupled thermo-mechanical eqns.

https://comet-fenics.readthedocs.io/en/latest/demo/thermoelasticity/thermoelasticity_transient.html#Variational-formulation-and-time-discretization

You can see within the time loop, it has defined
U = Function(V)
Why didn’t it work in my case?

5

Because if you inspect the variational formulation, you see that this is a linear problem, ie it is defined with a TrialFunction, and you are solving a==L, not F==0, which requires a Newton Method, see for instance
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html

6

Yeah makes sense. Thanks a lot.