Hello,
I am new to FEniCS and working towards solving nonlinear systems with the first-order generalized-alpha method (paper here). To get started, I’ve modified the dolfinx heat conduction tutorial (found here to solve with fem.petsc.NonlinearProblem() with generalized alpha. The solution looks reasonable at first glance, however, it does not produce the expected error convergence in the L^2-Norm.
Because the difference here is changing from linear to non-linear, I suspect there is a problem with how the NonlinearProblem is set up or how the nonlinear problem interprets the weak form. Both of these concern me because of the form that the Residual and Jacobian take for the generalized alpha method in terms of intermediate values of the solution variables. The code is below - thanks in advance!
import numpy
from dolfinx import mesh, fem, nls, plot
import ufl
from mpi4py import MPI
from petsc4py import PETSc
def run_problem(dt, rhoinf):
t = 0 # Start time
T = 2 # End time
num_steps = round(T/dt) # Number of time steps
# dt = (T-t)/num_steps # Time step size
alpha = 3
beta = 1.2
# generalized-alpha mehtod parameters
#--------------------------------------------------------------------------------
# rhoinf = 0.5 # rho infinity for the method, 1 corresponds to midpoint method
alpha_m = (1/2)*( (3-rhoinf)/(1+rhoinf))
alpha_f = 1/(1+rhoinf)
gamma = (1/2) + alpha_m - alpha_f
# domain setup
#--------------------------------------------------------------------------------
nx, ny = int(1/dt), int(1/dt)
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("CG", 1))
class exact_solution():
def __init__(self, alpha, beta, t):
self.alpha = alpha
self.beta = beta
self.t = t
def __call__(self, x):
return 1 + x[0]**2 + self.alpha * x[1]**2 + self.beta * self.t
class exact_du():
def __init__(self, alpha, beta, t):
self.alpha = alpha
self.beta = beta
self.t = t
def __call__(self, x):
return 0*(1 + x[0]**2 + self.alpha * x[1]**2) + self.beta * self.t
# Generalized-alpha helper functions
#--------------------------------------------------------------------------------
def predict_u(u_old):
return u_old
def predict_du(du_old, g):
return (g-1)/g * du_old
def correct_du(u_naf,u_old,du_old,am,af,g,dt):
return (1 - am/g)*du_old + (am/(g*dt*af))*(u_naf - u_old)
def update(old, mid, a):
return old + (1/a)*(mid - old)
# BCs
#--------------------------------------------------------------------------------
u_exact = exact_solution(alpha, beta, t)
du_exact = exact_du(alpha, beta, t)
u_D = fem.Function(V)
u_D.interpolate(u_exact)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))
# RHS function
#--------------------------------------------------------------------------------
f = fem.Constant(domain, beta - 2 - 2 * alpha)
#--------------------------------------------------------------------------------
# NONLINEAR SYSTEM SETUP WITH GEN ALPHA
#--------------------------------------------------------------------------------
# unknowns and intermediate variables
u_naf = fem.Function(V) # u at n + alpha_f
u_n = fem.Function(V)
du_n = fem.Function(V)
u_np1 = fem.Function(V)
du_np1 = fem.Function(V)
du_nam = fem.Function(V) # du/dt at n + alpha_m
# newmark formula
u_np1.x.array[:] = u_n.x.array + dt*du_n.x.array + gamma*dt*(du_np1.x.array - du_n.x.array)
du_nam.x.array[:] = correct_du(u_naf.x.array, u_n.x.array, du_n.x.array, alpha_m, alpha_f, gamma, dt)
# test function
v = ufl.TestFunction(V)
Res = du_nam*v*ufl.dx + ufl.dot(ufl.grad(u_naf), ufl.grad(v))*ufl.dx - f*v*ufl.dx
J = ufl.derivative(Res, u_naf, du_nam)
# setup non-linear problem
problem = fem.petsc.NonlinearProblem(Res,u_naf,[bc])
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = False
# initial conditions
u_n.interpolate(u_exact)
du_n.interpolate(du_exact)
# time looping
#--------------------------------------------------------------------------------
for n in range(num_steps):
# Update Diriclet boundary condition
u_exact.t+=alpha_f*dt # t at n+alpha_f
u_D.interpolate(u_exact)
# Solve NON-linear problem
#--------------------------------------------------------------------------------
# predict
u_np1.x.array[:] = predict_u(u_n.x.array)
du_np1.x.array[:] = predict_du(du_n.x.array,gamma)
# set intermediate values
u_naf.x.array[:] = u_n.x.array + alpha_f*( u_np1.x.array - u_n.x.array)
du_nam.x.array[:] = du_n.x.array + alpha_m*(du_np1.x.array - du_n.x.array)
# solver call to update u at n+alpha_f
r = solver.solve(u_naf)
# correct du/dt at n+alpha_m after Newton step
du_nam.x.array[:] = correct_du(u_naf.x.array, u_n.x.array, du_n.x.array, alpha_m, alpha_f, gamma, dt)
# update n+1
u_np1.x.array[:] = update( u_np1.x.array, u_naf.x.array, alpha_f)
du_np1.x.array[:] = update(du_np1.x.array, du_nam.x.array, alpha_m)
# update for next iteration
u_n.x.array[:] = u_np1.x.array
du_n.x.array[:] = du_np1.x.array
u_naf.x.scatter_forward()
u_exact.t+=(-alpha_f*dt + dt)
# error calcs
V2 = fem.FunctionSpace(domain, ("CG", 3))
uex = fem.Function(V2)
uex.interpolate(u_exact)
L2_error = fem.form(ufl.inner(u_n - uex, u_n - uex) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = numpy.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))
dof = len(u_naf.x.array)
return error_L2, u_naf, dof
import numpy
# looping parameters & variables
#--------------------------------------------------------------------------------
rhos = numpy.array([0, 0.5, 1])
dxs = numpy.array([0.2, 0.1, 0.05, 0.025])
L2_Errors = numpy.empty([3,4])
dofs = numpy.empty([3,4])
# run convergece study
#--------------------------------------------------------------------------------
for outer in range(3):
for inner in range(4):
# indexing rho infinity and time step/mesh size
rhoinf = rhos[outer]
dt = dxs[inner]
# run problem for convergence study
error_L2,u_naf,dof = run_problem(dt, rhoinf)
# save off L2-error and DOFs
L2_Errors[outer][inner] = error_L2
dofs[outer][inner] = dof # redundant
L2_Errors