Thank you for your response.
I am working on a time-dependent mixed formulation and used the solution here to formulate and solve the problem in the linear case.
Adding a nonlinearity to my formulation results in the Newton method not converging. In the code below I have labelled this as Method1.
My goal is to change the relaxation_parameter in each newton iteration.
Method2 follows the example you have given. Unfortunately I get
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
in each time step and I can’t explain why. It seem like the method update_solution doesn’t get called, when I run custom_solver.solve().
Another approach is Method3. Here I get the error
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
This post sasys the problem is resulting from the wrong use of split(). I tried both vuh.split() and split(vuh) but it doesnt make a difference
P1 = VectorElement('Lagrange', cell=triangle, degree=degree)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
tol = 1E-14
def clamped_boundary(x, on_boundary): # beam is only fixed on the left end
return on_boundary and near(x[0], 0.0, tol) #x[0] < tol
u_D1 = Constant((0, 0))
bc1 = DirichletBC(V.sub(0), u_D1, clamped_boundary)
bc2 = DirichletBC(V.sub(1), u_D1, clamped_boundary)
def force_on_rhs_boundary(x, on_boundary):
return on_boundary and near(x[0], 60.0, tol)
class topbot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and 0.0 + tol < x[0] < 60.0 - tol
kappa = 0.5
tt = np.linspace(0, 10, 100)
S = 1.0
SW = 0
fct = str(S) + '-(' + str(S) + '-' + str(SW) + ')*exp(-kappa*t)'
u_D22 = Expression((fct, '0'), degree=1, kappa=kappa, t=0)
u_D21 = Constant((0, 0))
bc3 = DirichletBC(V.sub(0), u_D21, force_on_rhs_boundary)
bc4 = DirichletBC(V.sub(1), u_D22, force_on_rhs_boundary)
# Mark boundaries as "ds"
boundaries = MeshFunction("size_t", mesh, dim=1)
boundaries.set_all(0)
#Mark boundaries where u*n=0
#tlt.mark(boundaries, 1)
topbot = topbot()
topbot.mark(boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
#Combine Dirichlet bcs
bcs = [bc3, bc4]
# Define strain and stress
def epsilon(z): # strain
return 0.5*(nabla_grad(z) + nabla_grad(z).T)
# return sym(nabla_grad(z))
def sigma(z): # stress
return lambda_*nabla_div(z) * Identity(d) + 2 * mu * epsilon(z)
# Define test functions
vu = Function(V)
v, u = split(vu)
vnun = Function(V)
v_n, u_n = split(vnun)
(phi1, phi2) = TestFunctions(V)
d = v.geometric_dimension()
# Define source terms
f_1 = Constant((0, 0))
f_2 = Constant((0, 0))
# Define variational problem
gamma = Constant(gamma_nitsche)
h = CellDiameter(mesh)
n = FacetNormal(mesh)
def ppos(x):
return (x+abs(x))/2
# https://github.com/MiroK/fenics-nitsche/blob/master/poisson/poisson_circle_dirichlet.py
# https://pastebin.com/3zPiJd9m # Bsp von 2018
# Weak formulation
F = (dot(v-v_n, phi2) * dx
- dt * inner(sigma(u), epsilon(phi2)) * dx
- dt * dot(v, phi1) * dx
+ dot(u-u_n, phi1) * dx)
# Add penalty term
F += gamma/h**2 * dot(ppos(dot(u, n)), dot(phi2, n)) * ds(1)
# Add Nitsche term:
F -= dot(dot(n, dot(sigma(u), n)), dot(phi2, n)) * ds(1)
class Problem(NonlinearProblem):
def __init__(self, a, L):
self.a = a
self.L = L
NonlinearProblem.__init__(self)
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
vuh = Function(V)
du = TrialFunction(V)
# Write custom Newton Solver:
class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())
def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)
PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps")
self.linear_solver().set_from_options()
def update_solution(self, x, dx, relaxation_parameter, nonlinear_problem,
iteration):
tau = 1.0
theta = min(sqrt(2.0*tau/norm(dx, norm_type="l2", mesh=V.mesh())), 1.0)
print(theta)
#info(f"Newton damping parameter: {theta:.3e}")
x.axpy(-theta, dx)
# Time stepping
for n in range(num_steps):
# Update current time
t += dt
print('t= ', t)
# Update bc
u_D22.t = t
J = derivative(F, vuh, du)
# Compute derivative of F for Newton method
R = action(F, vuh)
DR = derivative(R, vuh)
# Method1:
#problem = NonlinearVariationalProblem(R, vuh, bcs, DR)
#solver = NonlinearVariationalSolver(problem)
# Method2:
problem = Problem(DR, F)
custom_solver = CustomSolver()
custom_solver.solve(problem, vuh.vector())
# Method3:
"""
a_tol, r_tol = 1e-7, 1e-10
bcs_du = []
for bc in bcs:
bcs_du.append(bc)
print(bcs_du)
U_inc = Function(V)
nIter = 0
eps = 1
while eps > 1e-10 and nIter < 10: # Newton iterations
nIter += 1
A, b = assemble_system(J, -F, bcs_du)
solve(A, U_inc.vector(), b) # Determine step direction
eps = np.linalg.norm(U_inc.vector().array(), ord=2)
a = assemble(F)
for bc in bcs_du:
bc.apply(a)
#print 'b.norm("l2") = ', b.norm('l2'), 'np.linalg.norm(a.array(), ord = 2) = ', np.linalg.norm(a.array(), ord = 2)
#fnorm = b.norm('l2')
lmbda = 1.0 # step size, initially 1
vuh.vector()[:] += lmbda*U_inc.vector() # New u vector
print( {0:2d} {1:3.2E} {2:5e}'.format(nIter, eps, fnorm))
"""
# Split solutions
vh, uh = vuh.split(True)
# Update previous solution
v_n = vh
u_n = uh