Hi everyone,
I would like to use Newton’s method to solve a nonlinear problem in a mixed formulation, but I encountered the following error and would like to ask what might be missing in my code.
As far as I know, in DG methods, the mixed formulation is usually transformed into a primal formulation before solving. Therefore, I am unsure whether Newton’s method can be directly applied to the mixed formulation. Apologies, as I am not very familiar with this aspect.
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
import math
p = 4
pdeg = 1
size = 10
def vecnorm(vec):
return pow(inner(vec,vec),0.5)
def ConvergenceFunction(iteration, v, v0, abs_tol, rel_tol, convergence):
absolute = v.norm("l2")
if (iteration == 0):
absolute0 = absolute
else:
absolute0 = v0
relative = absolute/absolute0
if absolute < abs_tol or relative < rel_tol:
convergence = True
return convergence, absolute, relative
mesh = RectangleMesh(Point(0, 0), Point(1, 1), size, size)
uex = Expression(" sin(pi*x[0])*sin(pi*x[1])",degree=5,pi=math.pi)
Lex = Expression((" pi*sin(pi*x[1])*cos(pi*x[0]) "," pi*sin(pi*x[0])*cos(pi*x[1]) "),degree=5,pi=math.pi)
Coe = Expression("pow((pow(pi,2)*pow(sin(pi*x[0]),2)*pow(cos(pi*x[1]),2) + pow(pi,2)*pow(sin(pi*x[1]),2)*pow(cos(pi*x[0]),2)),1.0)", degree=5,pi=math.pi)
Aex = Expression((" Coe*pi*sin(pi*x[1])*cos(pi*x[0]) "," Coe*pi*sin(pi*x[0])*cos(pi*x[1]) "),degree=5,pi=math.pi,Coe=Coe)
f = Expression(" (0.5*pow(pi,4.0)*(-cos(pi*(2*x[0] - 2*x[1])) - cos(pi*(2*x[0] + 2*x[1])) + 2) - 2.0*pow(pi,4)*pow(cos(pi*x[0]),2)*cos(2*pi*x[1]) - 2.0*pow(pi,4)*cos(2*pi*x[0])*pow(cos(pi*x[1]),2))*sin(pi*x[0])*sin(pi*x[1])", degree=5,pi=math.pi)
P1 = VectorElement("DG", triangle, pdeg)
P2 = VectorElement("DG", triangle, pdeg)
P3 = FiniteElement("DG", triangle, pdeg)
P1P2P3 = MixedElement([P1, P2, P3])
W = FunctionSpace(mesh, P1P2P3)
w = Function(W)
(L,A,u) = w.split()
(X,Y,z) = TestFunctions(W)
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2
V1 = FunctionSpace(mesh, P1)
V2 = FunctionSpace(mesh, P2)
V3 = FunctionSpace(mesh, P3)
uD = uex
uD = project(uD, V3)
bc = DirichletBC(W.sub(2), uD, "on_boundary")
# bc = DirichletBC(W, uD, "on_boundary")
# Primal formulation and iter
delta = Constant(1e-10)
alpha = Constant(10.0)
u = project(uex, V3)
L = project(Lex, V1)
# u = un
# L = Ln
# A = project(Aex, V2)
Nonlinear1 = (delta + pow((vecnorm(L)),p-2))*L
Nonlinear2 = (delta + pow((vecnorm((1/h_avg)*jump(u,n))),p-2))*(1/h_avg)*jump(u,n)
Nonlinear3 = (delta + pow((vecnorm((1/h)*(u-uD)*n)),p-2))*(1/h)*(u-uD)
a = inner(Nonlinear1 , Y)*dx + alpha*inner(Nonlinear2,jump(z,n))*dS + alpha*Nonlinear3*z*ds
b1 = - inner(A,Y)*dx + inner(A,grad(z))*dx - inner(avg(A),jump(z,n))*dS - inner(A,n)*z*ds
b2 = - inner(X,L)*dx + inner(X,grad(u))*dx - inner(avg(X),jump(u,n))*dS - inner(X,n)*u*ds
F = f*z*dx
G = -uD*inner(X,n)*ds
Lhs = a + b1 + b2
Rhs = F + G
form = Lhs-Rhs
dw = TrialFunction(W)
dw_ = Function(W)
w_ = Function(W)
solve(form == 0, w, bc, J = derivative(form, u, du),
solver_parameters={"newton_solver":{"relative_tolerance":1e-6, "absolute_tolerance":1e-10, "convergence_criterion":"incremental"}})
*** Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /tmp/dolfin-src/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: 1c52e837eb54cc34627f90bde254be4aa8a2ae17