Occurence of error "Unable to successfully call PETSc function 'MatSetValuesLocal"

General
1

Dear. Researchers
I aim at solving two coupled time-dependent problems via the following code. upon execution of the code I have received the following error:

rror: 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-pkgs_1739259282872/work/dolfin/la/PETScMatrix.cpp.
*** Process: 0

I am not sure if this error stems from implementation or it possesses an external cause with PETSC. I would be thankful to assess the issue and give me your feedback.

import random
from dolfin import *
import ufl
from ufl.algorithms.formtransformations import compute_form_arities
from ufl.algorithms import expand_derivatives
import numpy as np
import matplotlib.pyplot as plt


dt     = 1e-1 # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson


mesh = RectangleMesh(Point(0,0), Point(200,200), 100,100)


P1 = FiniteElement("Lagrange", triangle, 1)
ME = FunctionSpace(mesh,MixedElement(P1, P1, P1))
mesh_file = File("Mesh1D.xml")
mesh_file << mesh
plt.figure()
plot(mesh, title = "domain")



# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step
du    = TrialFunction(ME)


v1 ,v2, v3 = TestFunctions(ME)
u1, u2, u3 = split(u)
r1, s1, t1 = split(u0)

du1, du2, du3 = split(du)

class K(UserExpression):
    
      
      
      def eval(self, values, x):
          c0 = 0.5
          eps = 0.01
          eps_eta = 0.1
          Psi = 1.5
          
          values[0] = cos(x[0])*cos(x[1])
          
          
         
          values[1] = 0.0  
          values[2] = cos(x[0]) * cos(x[0])
         
      def value_shape(self):
          return(3, )
   
        
order_value = K(mesh)



u = interpolate(order_value, ME)


u0.assign(u)

#define parameters

c_alpha = Constant(0.3) 
c_beta = Constant(0.7)
ro_s = sqrt(2)
M = Constant(5)
L = Constant(5)
kappa_c = Constant(3)
kappa_eta = Constant(3)
w = Constant(1)
alpha = Constant(5)



u1 = variable(u1)

f_a = (ro_s**2)*(u1 - c_alpha)**2

f_b = (ro_s**2)*(c_beta - u1)**2

fder_a = diff(f_a , u1)
fder_b = diff(f_b, u1)
# define swtiching function

u3 = variable(u3)


hfunction = (u3**3*(6*u3**2 - 15*u3 + 10))



g = ((u3**2)*(1-u3)**2)

                                   
fchem = (f_a * (1 - hfunction) + f_b * hfunction + w*g)
      
fchem_c = diff(fchem,u1)
fchem_eta1 = diff(fchem,u3)


Form_cons= u1*v1*dx - r1*v1*dx + dt*M*inner(grad(u2), grad(v1))*dx



Form_mu = u2*v2*dx - fchem_c*v2*dx - kappa_c*inner(grad(u1), grad(v2))*dx 


Form_eta1 = (u3*v3*dx - t1*v3*dx + dt*L*fchem_eta1*v3*dx + 
            L*kappa_eta*dt*inner(grad(u3), grad(v3))*dx)
            
Form = Form_cons + Form_mu + Form_eta1


Gain = derivative(Form, u, du)


t = 0.0
T = 1000*dt
while (t < T):
    t += dt
    solve(Form == 0, u , J = Gain)
    u0.assign(u)
2

You have defined u twice.
Both here

The latter one should be u.interpolate(order_value)
Otherwise you differentiate F with respect to a variable that is not in your variational form.

3

Dear. Dr. Dokken
Thanks for your prompt response. After applying modification, the code worked but the convergence behavior was not achieved. To remedy the issue, I changed the nonlinear solver to SNES but not effective. Since I adopted all the parameter from a benchmark I expect normal solution with good convergence behavior. I am so grateful if you assess the updated code and give me your feedback.

import random
from dolfin import *
import ufl
from ufl.algorithms.formtransformations import compute_form_arities
from ufl.algorithms import expand_derivatives
import numpy as np
import matplotlib.pyplot as plt


dt     = 1e-1 # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson


mesh = RectangleMesh(Point(0,0), Point(200,200), 100,100)



P1 = FiniteElement("Lagrange", triangle, 1)
ME = FunctionSpace(mesh,MixedElement(P1, P1, P1))
mesh_file = File("Mesh1D.xml")
mesh_file << mesh
plt.figure()
plot(mesh, title = "domain")



# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step
du    = TrialFunction(ME)


v1 ,v2, v3 = TestFunctions(ME)
u1, u2, u3 = split(u)
r1, s1, t1 = split(u0)

du1, du2, du3 = split(du)

class K(UserExpression):
    
      
      
      def eval(self, values, x):
          c0 = 0.5
          eps = 0.01
          eps_eta = 0.1
          Psi = 1.5
          
          values[0] = cos(x[0]) * cos(x[1])
          values[1] = 0.0  
          values[2] = cos(x[0])*cos(x[1])
         
      def value_shape(self):
          return(3, )
   
        
order_value = K(mesh)



u.interpolate(order_value)



file1 = File("xglobal_ini.pvd")
file2 = File("chemical potential.pvd")
file3 = File("eta1_ini.pvd")


file1 << u.split(deepcopy= True)[0]

file2 << u.split(deepcopy= True)[1]
file3 << u.split(deepcopy= True)[2]



u0.assign(u)

#define parameters

c_alpha = 0.3
c_beta = 0.7
ro_s = sqrt(2)
M = Constant(5)
L = Constant(5)
kappa_c = Constant(3)
kappa_eta = Constant(3)
w = Constant(1)
alpha = Constant(5)



u1 = variable(u1)

f_a = (ro_s**2)*(u1 - c_alpha)**2

f_b = (ro_s**2)*(c_beta - u1)**2
fder_a = diff(f_a , u1)
fder_b = diff(f_b, u1)


# define swtiching function

u3 = variable(u3)


hfunction = (u3**3*(6*u3**2 - 15*u3 + 10))

g = ((u3**2)*(1-u3)**2)



                                   
fchem = (f_a * (1 - hfunction) + f_b * hfunction + w*g)
      
fchem_c = diff(fchem,u1)
fchem_eta1 = diff(fchem,u3)




Form_cons = u1*v1*dx - r1*v1*dx + dt*M*inner(grad(u2), grad(v1))*dx



Form_mu = u2*v2*dx - fchem_c*v2*dx - kappa_c*inner(grad(u1), grad(v2))*dx 


Formt = (u3*v3*dx - t1*v3*dx +
         dt*L*fchem_eta1*v3*dx + 
            L*kappa_eta*dt*inner(grad(u3), grad(v3))*dx)


            
Form = Form_cons + Form_mu + Formt


Gain = derivative(Form, u, du)



problem = NonlinearVariationalProblem(Form, u, J = Gain)
solver = NonlinearVariationalSolver(problem)

solver.parameters['nonlinear_solver'] = 'snes'
solver.parameters['snes_solver']['linear_solver'] = 'mumps'

solver.parameters['snes_solver']['report'] = True
solver.parameters['snes_solver']['line_search'] = 'bt'

solver.parameters['snes_solver']['error_on_nonconvergence'] = False
solver.parameters['snes_solver']['absolute_tolerance'] = 1e-6


t = 0.0
T = 1000*dt
while (t < T):
    t += dt
    solver.solve()
    u0.assign(u)