PETSc error 63 for a simple nonlinear 1D simulation

Hello!

I got this error for my code.

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

After look up for some posts in this discourse on the PETSc error 63, I am still unable to solve this issue for a really simple demo. My code is really simple and it wants to solve two coupled nonlinear PDEs.

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)

W_elem = MixedElement([Element1, Element2])
ME = FunctionSpace(mesh, W_elem)

# Define functions
u_old = Function(ME)
u_new = Function(ME)

Qx_old, Qy_old= split(u_old)
Qx_new, Qy_new = split(u_old)
Qx_mid = 0.5 * Qx_old + 0.5 * Qx_new
Qy_mid = 0.5 * Qy_old + 0.5 * Qy_new
tf = TestFunction(ME)
q, v = split(tf)
dt = 0.01

# Neumann Boundary Condition
def boundary(x, on_boundary):
    return on_boundary


bc_N = DirichletBC(ME.sub(1), Constant(0), boundary)

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)


    def value_shape(self):
        return (2,)

# Set initial conditions to old function
u_init = InitialConditions()
u_old.interpolate(u_init)

# Weak statement of the equations
F1 = (Qx_new - Qx_old) / dt * q * dx - Qy_mid * Qx_mid.dx(0) * q * dx
F0 = (Qy_new - Qy_old) / dt * v * dx - Qx_mid ** 2 * v * dx

F = F0 + F1

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()

How can I modify this code in order to run it successfully? I thought the solutions might be just trivial, but even in that case, can I somehow rewrite this code to solve it?

Any help and advice is greatly appreciated!

Qx_new, Qy_new = split(u_old)

If you change this to

Qx_new, Qy_new = split(u_new)

do you actually get a non-zero Jacobian?

2 Likes

Hi nate!

Thank you so much. I just realized that I made such a silly mistake haha. I was so confused yesterday and try a lot things. Fixing this typo makes the code running. Thanks again! Have a good day!