Please note that you are not enforcing your force term correctly.
This term is not set in your variational form.
This is how I would write your code:
class ManufacturedSolution:
def __init__(self):
self.t = 0
def eval(self, x):
values = np.sin(2*np.pi*self.t) * np.sin(2*np.pi*x[0]) * np.cos(2*np.pi*x[1])
return values
import numpy as np
from dolfinx import mesh, fem, nls, plot, cpp, io
import ufl
from mpi4py import MPI
from petsc4py import PETSc
### mesh and problem parameters
#--------------------------------------------------------------------------------
dt = 0.05
nstep = int(2/dt)
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0., 0.), (1., 1.)),n=(20, 20), cell_type=mesh.CellType.triangle)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
n = ufl.FacetNormal(domain)
x = ufl.SpatialCoordinate(domain)
filename = 'Mfg.xdmf'
with io.XDMFFile(domain.comm, filename, 'w') as xdmf:
xdmf.write_mesh(domain)
#--------------------------------------------------------------------------------
# weak form setup
#--------------------------------------------------------------------------------
# functions
W = fem.FunctionSpace(domain, ("CG",1))
u_n = fem.Function(W)
u_np1 = fem.Function(W)
v = ufl.TestFunction(W)
dU = (u_np1 - u_n)/dt
t_nh = fem.Constant(domain, dt/2)
force = 2*ufl.pi*ufl.cos(2*ufl.pi*t_nh)*ufl.sin(2*ufl.pi*x[0])*ufl.cos(2*ufl.pi*x[1])
Res = ufl.inner( dU, v)*ufl.dx - ufl.inner(force, v)*ufl.dx
#--------------------------------------------------------------------------------
# BCs
#--------------------------------------------------------------------------------
solFunc = ManufacturedSolution()
solFunc.t = 0
bcFunc = fem.Function(W)
bcFunc.name = 'Analytic'
bcFunc.interpolate(solFunc.eval)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(W, fdim, boundary_facets)
bc = fem.dirichletbc(bcFunc, boundary_dofs)
bcs = [bc]
#--------------------------------------------------------------------------------
# Nonlinear problem setup
#--------------------------------------------------------------------------------
J = ufl.derivative(Res, u_np1)
residual = fem.form(Res)
jacobian = fem.form(J)
# time dependent structures
A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(domain.comm)
opts = PETSc.Options()
prefix = f"solver_{id(solver)}"
solver.setOptionsPrefix(prefix)
option_prefix = solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
# opts[f"{option_prefix}ksp_monitor"] = None # uncomment to view linear solver output
opts[f"{option_prefix}pc_type"] = "lu"
solver.setFromOptions()
solver.setOperators(A)
DeltaA = fem.Function(W)
max_iterations = 5
# initial condition
#--------------------------------------------------------------------------------
ICFunc = fem.Function(W)
ICFunc.interpolate(solFunc.eval)
u_n.x.array[:] = ICFunc.x.array
u_np1.x.array[:] = ICFunc.x.array
with io.XDMFFile(domain.comm, filename, 'a') as xdmf:
xdmf.write_function(u_np1, 0)
xdmf.write_function(bcFunc,0)
#--------------------------------------------------------------------------------
# Time Loop
#--------------------------------------------------------------------------------
t_n = 0
bc = fem.dirichletbc(bcFunc, boundary_dofs)
bcs = [bc]
for niter in range(nstep):
# update time dependent force and BCs
t_nh.value = t_n + dt/2
solFunc.t = t_nh.value+dt
bcFunc.interpolate(solFunc.eval)
# Solve NON-linear problem
#--------------------------------------------------------------------------------
# begin newton loop
i = 0
u_np1.x.array[:] = u_n.x.array
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, jacobian, bcs = bcs)
A.assemble()
fem.petsc.assemble_vector(L, residual)
L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
# Scale residual by -1
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [jacobian], [bcs], x0=[u_np1.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, [bc], u_np1.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linearized problem
solver.solve(L, DeltaA.vector)
DeltaA.x.scatter_forward()
# Update du_{i+1} = du_i + delta x_i
u_np1.x.array[:] += DeltaA.x.array
i += 1
# Compute norm of update
correction_norm = DeltaA.vector.norm(0)
print(f"Time step {niter}, Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
# output
with io.XDMFFile(domain.comm, filename, 'a') as xdmf:
xdmf.write_function(u_np1, t_n+dt)
xdmf.write_function(bcFunc, t_n+dt)
# update for next time step
u_n.x.array[:] = u_np1.x.array
t_n+=dt
It gives the correct bc, but the wrong force.
To debug such things (including the time stepiing), I would strongly suggest first solving the linear problem:
class ManufacturedSolution:
def __init__(self):
self.t = 0
def eval(self, x):
values = np.sin(2*np.pi*self.t) * np.sin(2*np.pi*x[0]) * np.cos(2*np.pi*x[1])
return values
import numpy as np
from dolfinx import mesh, fem, io
import ufl
from mpi4py import MPI
### mesh and problem parameters
#--------------------------------------------------------------------------------
dt = 0.1
nstep = int(2/dt)
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0., 0.), (1., 1.)),n=(20, 20), cell_type=mesh.CellType.triangle)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
n = ufl.FacetNormal(domain)
x = ufl.SpatialCoordinate(domain)
filename = 'Mfg.xdmf'
with io.XDMFFile(domain.comm, filename, 'w') as xdmf:
xdmf.write_mesh(domain)
#--------------------------------------------------------------------------------
# weak form setup
#--------------------------------------------------------------------------------
# functions
W = fem.FunctionSpace(domain, ("CG",1))
u_n = fem.Function(W)
u_np1 = fem.Function(W)
v = ufl.TestFunction(W)
dU = (u_np1 - u_n)
t_nh = fem.Constant(domain, dt/2)
force = 2*ufl.pi*ufl.cos(2*ufl.pi*t_nh)*ufl.sin(2*ufl.pi*x[0])*ufl.cos(2*ufl.pi*x[1])
Res = ufl.inner( dU, v)*ufl.dx - dt*ufl.inner(force, v)*ufl.dx
uh = ufl.TrialFunction(W)
F = ufl.replace(Res, {u_np1:uh})
a, L = ufl.system(F)
#--------------------------------------------------------------------------------
# BCs
#--------------------------------------------------------------------------------
solFunc = ManufacturedSolution()
solFunc.t = 0
bcFunc = fem.Function(W)
bcFunc.name = 'Analytic'
bcFunc.interpolate(solFunc.eval)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(W, fdim, boundary_facets)
bc = fem.dirichletbc(bcFunc, boundary_dofs)
bcs = [bc]
lin_solver = fem.petsc.LinearProblem(a, L, u=u_np1, bcs=bcs)
# initial condition
#--------------------------------------------------------------------------------
ICFunc = fem.Function(W)
ICFunc.interpolate(solFunc.eval)
u_n.x.array[:] = ICFunc.x.array
u_np1.x.array[:] = ICFunc.x.array
with io.XDMFFile(domain.comm, filename, 'a') as xdmf:
xdmf.write_function(u_np1, 0)
xdmf.write_function(bcFunc,0)
#--------------------------------------------------------------------------------
# Time Loop
#--------------------------------------------------------------------------------
t_n = 0
bc = fem.dirichletbc(bcFunc, boundary_dofs)
bcs = [bc]
for niter in range(nstep):
# update time dependent force and BCs
t_nh.value = t_n + dt/2
solFunc.t = t_nh.value+dt/2
bcFunc.interpolate(solFunc.eval)
lin_solver.solve()
u_n.x.array[:] = u_np1.x.array
t_n += dt
with io.XDMFFile(domain.comm, filename, 'a') as xdmf:
xdmf.write_function(u_np1, t_n)
xdmf.write_function(bcFunc,t_n)
which makes it easy to fix the non-linear version
class ManufacturedSolution:
def __init__(self):
self.t = 0
def eval(self, x):
values = np.sin(2*np.pi*self.t) * np.sin(2*np.pi*x[0]) * np.cos(2*np.pi*x[1])
return values
import numpy as np
from dolfinx import mesh, fem, nls, plot, cpp, io
import ufl
from mpi4py import MPI
from petsc4py import PETSc
### mesh and problem parameters
#--------------------------------------------------------------------------------
dt = 0.1
nstep = int(2/dt)
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0., 0.), (1., 1.)),n=(20, 20), cell_type=mesh.CellType.triangle)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
n = ufl.FacetNormal(domain)
x = ufl.SpatialCoordinate(domain)
filename = 'Mfg.xdmf'
with io.XDMFFile(domain.comm, filename, 'w') as xdmf:
xdmf.write_mesh(domain)
#--------------------------------------------------------------------------------
# weak form setup
#--------------------------------------------------------------------------------
# functions
W = fem.FunctionSpace(domain, ("CG",1))
u_n = fem.Function(W)
u_np1 = fem.Function(W)
v = ufl.TestFunction(W)
dU = (u_np1 - u_n)/dt
t_nh = fem.Constant(domain, dt/2)
force = 2*ufl.pi*ufl.cos(2*ufl.pi*t_nh)*ufl.sin(2*ufl.pi*x[0])*ufl.cos(2*ufl.pi*x[1])
Res = ufl.inner( dU, v)*ufl.dx - ufl.inner(force, v)*ufl.dx
#--------------------------------------------------------------------------------
# BCs
#--------------------------------------------------------------------------------
solFunc = ManufacturedSolution()
solFunc.t = 0
bcFunc = fem.Function(W)
bcFunc.name = 'Analytic'
bcFunc.interpolate(solFunc.eval)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(W, fdim, boundary_facets)
bc = fem.dirichletbc(bcFunc, boundary_dofs)
bcs = [bc]
#--------------------------------------------------------------------------------
# Nonlinear problem setup
#--------------------------------------------------------------------------------
J = ufl.derivative(Res, u_np1)
residual = fem.form(Res)
jacobian = fem.form(J)
# time dependent structures
A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(domain.comm)
opts = PETSc.Options()
prefix = f"solver_{id(solver)}"
solver.setOptionsPrefix(prefix)
option_prefix = solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
# opts[f"{option_prefix}ksp_monitor"] = None # uncomment to view linear solver output
opts[f"{option_prefix}pc_type"] = "lu"
solver.setFromOptions()
solver.setOperators(A)
DeltaA = fem.Function(W)
max_iterations = 5
# initial condition
#--------------------------------------------------------------------------------
ICFunc = fem.Function(W)
ICFunc.interpolate(solFunc.eval)
u_n.x.array[:] = ICFunc.x.array
u_np1.x.array[:] = ICFunc.x.array
with io.XDMFFile(domain.comm, filename, 'a') as xdmf:
xdmf.write_function(u_np1, 0)
xdmf.write_function(bcFunc,0)
#--------------------------------------------------------------------------------
# Time Loop
#--------------------------------------------------------------------------------
t_n = 0
bc = fem.dirichletbc(bcFunc, boundary_dofs)
bcs = [bc]
for niter in range(nstep):
# update time dependent force and BCs
t_nh.value = t_n + dt/2
solFunc.t = t_n + dt
bcFunc.interpolate(solFunc.eval)
# Solve NON-linear problem
#--------------------------------------------------------------------------------
# begin newton loop
i = 0
u_np1.x.array[:] = u_n.x.array
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, jacobian, bcs = bcs)
A.assemble()
fem.petsc.assemble_vector(L, residual)
L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
# Scale residual by -1
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [jacobian], [bcs], x0=[u_np1.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, [bc], u_np1.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linearized problem
solver.solve(L, DeltaA.vector)
DeltaA.x.scatter_forward()
# Update du_{i+1} = du_i + delta x_i
u_np1.x.array[:] += DeltaA.x.array
i += 1
# Compute norm of update
correction_norm = DeltaA.vector.norm(0)
print(f"Time step {niter}, Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
# output
with io.XDMFFile(domain.comm, filename, 'a') as xdmf:
xdmf.write_function(u_np1, t_n+dt)
xdmf.write_function(bcFunc, t_n+dt)
# update for next time step
u_n.x.array[:] = u_np1.x.array
t_n+=dt