Hi everyone. I meet an issue with my code while solving five time dependent PDEs. Any help and suggestions will be appreciated. I got the following error:
*** Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
I have already looked up for some posts with the same error. However, my situation is different and I am still unable to solve this issue even though I tried some of them.
The code is pretty straight forward. I have a lot of derivatives since I have five variables. I just write the five equations in F0, F1, F2, F3, F4, and add them up to F, then use solve(F == 0, u_new, bc_N) to solve it
from dolfin import *
from tqdm import tqdm
import random
import numpy as np
mesh = IntervalMesh(10000, 0, 1)
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element2 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element3 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element4 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element5 = FiniteElement("CG", mesh.ufl_cell(), 1)
W_elem = MixedElement([Element1, Element2, Element3, Element4, Element5])
ME = FunctionSpace(mesh, W_elem)
# Define functions
u_old = Function(ME)
u_new = Function(ME)
Qxx_old, Qxy_old, Cxx_old, Cxy_old, Vy_old = split(u_old)
Qxx_new, Qxy_new, Cxx_new, Cxy_new, Vy_new = split(u_old)
Qxx_mid = 0.5 * Qxx_old + 0.5 * Qxx_new
Qxy_mid = 0.5 * Qxy_old + 0.5 * Qxy_new
Cxx_mid = 0.5 * Cxx_old + 0.5 * Cxx_new
Cxy_mid = 0.5 * Cxy_old + 0.5 * Cxy_new
Vy_mid = 0.5 * Vy_old + 0.5 * Vy_new
tf = TestFunction(ME)
q, v, w, n, m = split(tf)
# Define parameters
lam = 1.5
ga_Q = 1
k_ga_Q = 1
kappa_ga_Q = 1
Ga = 1
eta = 1
alpha = -500
G = 1
beta = 1
kappa = 1
tau = 1
k_ga_p = 1
dt = 0.01
# Neumann Boundary Condition
def boundary(x, on_boundary):
return on_boundary
bc_N = DirichletBC(ME.sub(4), Constant(0), boundary)
# Initial Condition
class InitialConditions(UserExpression):
def eval(self, values, x):
values[0] = - 0.5 + 0.01 * random.gauss(0, 1)
values[1] = 0.01 * random.gauss(0, 1)
values[2] = 0.01 * random.gauss(0, 1)
values[3] = 0.01 * random.gauss(0, 1)
values[4] = 0
def value_shape(self):
return (5,)
# Set initial conditions to old function
u_init = InitialConditions()
u_old.interpolate(u_init)
# Weak statement of the equations
F0 = (Qxx_new - Qxx_old) / dt * q * dx - Qxy_mid * Vy_mid.dx(0) * q * dx - 1 / ga_Q * q * dx + \
4 * Qxx_mid ** 3 / ga_Q * q * dx + 4 * Qxy_mid ** 2 * Qxx_mid / ga_Q * q * dx + \
k_ga_Q * Qxx_mid.dx(0) * q.dx(0) * dx + 2 * kappa_ga_Q * Cxx_mid * q * dx
F1 = (Qxy_new - Qxy_old) / dt * v * dx - Qxx_mid * Vy_mid.dx(0) * v * dx - lam / 2 * Vy_mid.dx(0) * v * dx \
- 1 / ga_Q * v * dx + 4 * Qxy_mid ** 3 / ga_Q * v * dx + 4 * Qxx_mid ** 2 * Qxy_mid / ga_Q * v * dx + \
k_ga_Q * Qxy_mid.dx(0) * v.dx(0) * dx + 2 * kappa_ga_Q * Cxy_mid * v * dx
F2 = Ga * Vy_new * w * dx + eta * Vy_new.dx(0) * w.dx(0) * dx - alpha * Qxy_new.dx(0) * w * dx - \
G * Cxy_new.dx(0) * w * dx
F3 = (Cxx_new - Cxx_old) / dt * n * dx + Cxy_mid * Vy_mid.dx(0) * n * dx + 1 / tau * Cxx_mid * n * dx \
+ 2 * k_ga_p * Qxx_mid * n * dx
F4 = (Cxy_new - Cxy_old) / dt * m * dx - Cxx_mid * Vy_mid.dx(0) * m * dx + 1 / tau * Cxy_mid * m * dx \
- Vy_mid.dx(0) * m * dx + 2 * k_ga_p * Qxy_mid * m * dx
F = F0 + F1 + F2 + F3 + F4
t = 0.0
T = 1000 * dt
for t in tqdm(np.arange(0, 10.01, 0.01)):
solve(F == 0, u_new, bc_N)
u_old.vector()[:] = u_new.vector()
u_copy = Function(ME)
u_copy.vector()[:] = u_new.vector()