Time-dependent Dirichlet BCs are violated

[dg0609]

Hello,

I’m observing that my simple manufactured problem with a custom Newton solver has Dirichlet boundary conditions that are violated. The manufactured solution on the unit square is given by \phi = sin(2 \pi t) sin(2 \pi x) cos(2 \pi y) , and the PDE is simplified to \frac{\partial \phi}{\partial t} =2 \pi cos(2 \pi t) sin(2 \pi x) cos(2 \pi y) . The time-dependent Dirichlet BCs are defined from the manufactured solution.

The corresponding weak form residual is ufl.inner( dU, v)*ufl.dx - ufl.inner(force, v)*ufl.dx where 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]) and dU = (u_np1 - u_n)/dt .

This is what I’m seeing (note the zero values around the boundaries)

and this is what I expect to see (output from the actual function that is implemented to specify the Dirichlet BCs).

This is the MWE

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
import pyvista

### 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 = 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
for niter in range(nstep):

    # update time dependent force and BCs
    t_nh = t_n + dt/2
    solFunc.t = t_nh+dt
    bcFunc.interpolate(solFunc.eval)
    bc = fem.dirichletbc(bcFunc, boundary_dofs)
    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])

    # 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 = [bc])
        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], [[bc]], 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
        

There are many examples in the tutorials that I followed, so this is puzzling me. I welcome any feedback.

Thanks!

[dokken]

Please note that you are not enforcing your force term correctly.

@dg0609

for niter in range(nstep):

    # update time dependent force and BCs
    t_nh = t_n + dt/2
    solFunc.t = t_nh+dt
    bcFunc.interpolate(solFunc.eval)
    bc = fem.dirichletbc(bcFunc, boundary_dofs)
    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])

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

[ampdes]

I have a conceptual follow up question on this.

When solving a linear time dependent PDE problem, the computational graph is set and accessed to update the right hand side data at each time point via the functions used in defining the variational forms.

Similarly, I guess the function used in dirichletbc can be accessed to update the time varying dirichlet boundary conditons.

But what is

fem.petsc.assemble_matrix(MATRIX, FORM, bcs=bcs) doing?

For a linear problem, I guess it only needs to know the dirichlet boundary dofs, hence need not be assembled at each time point?

Thanks.

[dokken]

@ampdes

For a linear problem, I guess it only needs to know the dirichlet boundary dofs, hence need not be assembled at each time point?

Correct, we use lifting, ie, we assemble a matrix A where all Dirichlet dof rows and columns are set to the identity row/column. Then L is corrected with lifting, by computing L-= Ag Where g is equal to the Dirichlet bc value at those dofs, zero everywhere else.

Ive posted a reference for lifting at Difference between LUSolver and solve - #2 by dokken
where I Also highlight the benefit of lifting, ie symmetry is preserved.