Solving dynamic linear system without integrating the linear form

dolfinx
, ,
1

This question is strongly related to this other thread. So the main idea is:

I’m triying to couple DOLFINx via preCICE by myself, because the current adapter can solve 3D linear elastic problems at the moment.

In that threat, I solve the coupling problems I can’t apply loads over a DoF by directly manipulating the linear form, considering that FenicsX does not have a PointSource implementation I’m trying to use this approach.

I have prepared the following MWE that I test against analytical and CaluliX results accurately. It represents a long cylindrical beam (same as in the other thread) 150Lx6.5D, with a load of 10 in the direction of the X axis, on a DoF near the upper extreme, and clamped on its base.

import os
from typing import NamedTuple

# enforce 1 thread per core
nthreads = 1
os.environ["OMP_NUM_THREADS"] = str(nthreads)
os.environ["OPENBLAS_NUM_THREADS"] = str(nthreads)
os.environ["MKL_NUM_THREADS"] = str(nthreads)

import gmsh
import numpy as np
import ufl
from dolfinx import fem, io, mesh
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem.petsc import assemble_matrix, set_bc

comm = MPI.COMM_WORLD
CURRENT_FOLDER = os.path.dirname(__file__)
RESULTS_DIR = os.path.join(CURRENT_FOLDER, "results")
os.makedirs(RESULTS_DIR, exist_ok=True)
os.chdir(CURRENT_FOLDER)

# -------------------#
# Problem properties #
# -------------------#
class Material(NamedTuple):
    name: str
    E: PETSc.ScalarType
    nu: PETSc.ScalarType
    rho: PETSc.ScalarType


material = Material("Aluminum 6061", 68.9e9, 0.33, 2700)

# beam size
height = 150  # m
diameter = 6.5  # m

# loads
load = np.array([10.0, 0, 0])
# --------------------#


# ------- #
# MESHING #
# ------- #
FIXED_TAG = 1000
BEAM_SURFACE_TAG = 2000
TIP_SURFACE_TAG = 3000
VOLUME_TAG = 4000
ELEMENTS_ORDER = 2
VECTOR_SPACE_DEGREE = 2


def gmsh_tower(model: gmsh.model, name: str, h: float, d: float) -> gmsh.model:
    gmsh.clear()
    model.add(name)
    model.setCurrent(name)

    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.05)
    gmsh.option.setNumber("Mesh.MshFileVersion", 2)

    circle = model.occ.addDisk(0, 0, 0, d / 2, d / 2)
    surface = model.occ.addPlaneSurface([circle])
    model.occ.synchronize()

    extruded_geometry = model.occ.extrude([(2, circle)], 0, 0, h, numElements=[h], recombine=True)
    model.occ.synchronize()

    base_surf_pg = model.addPhysicalGroup(2, [surface], tag=FIXED_TAG, name="FIXED_SURFACE")

    subdivision = [h]
    extrusion = model.occ.extrude([(2, surface)], 0, 0, h, subdivision)
    model.occ.synchronize()

    volume = model.addPhysicalGroup(3, [extrusion[1][1]], tag=VOLUME_TAG, name="VOLUME")
    lateral_surf_group = model.addPhysicalGroup(
        2, [extrusion[2][1]], tag=BEAM_SURFACE_TAG, name="BEAM_SURFACE"
    )
    upper_surf_group = model.addPhysicalGroup(
        2, [extrusion[0][1]], tag=TIP_SURFACE_TAG, name="TIP_SURFACE"
    )

    model.mesh.generate(1)
    model.mesh.generate(2)
    model.mesh.generate(3)
    model.mesh.setOrder(ELEMENTS_ORDER)
    model.mesh.optimize()
    return model


def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, mode: str):
    msh, ct, ft = io.gmshio.model_to_mesh(model, comm, rank=0)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    return (msh, ct, ft)


gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model_name = "Tower"
model = gmsh.model()
model = gmsh_tower(model=model, name=model_name, h=height, d=diameter)
model.setCurrent(model_name)
domain, cell_markers, facet_markers = create_mesh(
    comm=MPI.COMM_SELF,
    model=model,
    name=model_name,
    mode="w",
)

# -------------- #
# Function Space #
# -------------- #
dim = domain.geometry.dim
V = fem.functionspace(domain, ("Lagrange", VECTOR_SPACE_DEGREE, (dim,)))

E = fem.Constant(domain, float(material.E))
nu = fem.Constant(domain, float(material.nu))
rho = fem.Constant(domain, float(material.rho))

# --------------------#
# Boundary conditions #
# --------------------#
# Fixed BC
fixed_facets = facet_markers.find(FIXED_TAG)
fixed_surface_dofs = fem.locate_dofs_topological(V, 2, fixed_facets)

u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, fixed_surface_dofs)]

# -------------------------#
# linear elastic equations #
# -------------------------#

lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)


def epsilon(u):
    return ufl.sym(ufl.grad(u))


def sigma(u):
    return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


def dof_points_in_vector_space(V, fn):
    bs = V.dofmap.bs

    fs_coords = V.tabulate_dof_coordinates()
    boundary_dofs = fem.locate_dofs_geometrical(V, fn)
    unrolled_dofs = np.empty(len(boundary_dofs) * bs, dtype=np.int32)

    for i, dof in enumerate(boundary_dofs):
        for b in range(bs):
            unrolled_dofs[i * bs + b] = dof * bs + b

    return unrolled_dofs.reshape(-1, bs), fs_coords[boundary_dofs]


def point(x):
    WP = (4.9490529e-02, -1.4846031e-01, 1.5000000e02)
    return np.isclose(x[0], WP[0]) & np.isclose(x[1], WP[1]) & np.isclose(x[2], WP[2])


point_dofs, point_dofs_coord = dof_points_in_vector_space(V, point)

du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
a = fem.form(ufl.inner(sigma(du), epsilon(u_)) * ufl.dx)

u = fem.Function(V, name="Displacement")
f = fem.Function(V, name="Load")

A = assemble_matrix(a, bcs=bcs)
A.assemble()

f.x.array[:] = 0
f.vector[point_dofs[0]] = load
f.vector.assemble()
set_bc(f.vector, bcs)


solver = PETSc.KSP().create(domain.comm)  # type: ignore
solver.setFromOptions()
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.CG)
pc = solver.getPC()
pc.setType(PETSc.PC.Type.GAMG)
solver.setTolerances(rtol=1e-5, atol=1e-11, max_it=300)

# Set a monitor, solve linear system, and display the solver
# configuration
solver.setConvergenceTest(lambda ksp, n, rnorm: print(f"Iteración {n}, Residuo {rnorm}"))
solver.solve(f.vector, u.vector)

# Scatter forward the solution vector to update ghost values
u.x.scatter_forward()

print("====================================================================")
print("                           DEFORMATION                              ")
print("====================================================================")
R = diameter / 2
L = height
Fx, Fy, Fz = load
Ix = np.pi / 4 * R**4
print(f"         Analytic: {(Fx * L**3) / (3 * material.E * Ix)}")
print(f"         Computed: {u.x.array[point_dofs]}")

So as this works fine, I tried to implement the FSI classic code tutorial perpendicular flap which was also implemented for the legacy version of fenics.

First I tried this code, In which I integrated the linear form over the coupling interface.

import numpy as np
import ufl
import os
from mpi4py import MPI

import dolfinx as dfx
import dolfinx.fem.petsc
from dolfinx.mesh import create_rectangle, CellType
from petsc4py import PETSc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, set_bc
from dolfinx.fem.petsc import LinearProblem
from turbine.solvers.interface import Adapter
from turbine.materials import Material
# --------- #
# CONSTANTS #
# --------- #

MPI_COMM = MPI.COMM_WORLD
CURRENT_FOLDER = os.path.dirname(__file__)
PARTICIPANT_CONFIG = os.path.join(CURRENT_FOLDER, "precice-config.json")
RESULTS_DIR = os.path.join(CURRENT_FOLDER, "results")
os.makedirs(RESULTS_DIR, exist_ok=True)
os.chdir(CURRENT_FOLDER)

WRITER = dfx.io.VTKFile(MPI_COMM, f"{RESULTS_DIR}/result.pvd", "w")

WIDTH, HEIGHT = 0.1, 1 
NX, NY = 4, 26


E = 4000000.0
NU = 0.3
RHO = 3000.0
material =  Material("Flap",E=E,nu=NU,rho=RHO)

BETA_ = 0.25
GAMMA_ = 0.5

# ------- #
# MESHING #
# ------- #

domain = create_rectangle(
    MPI_COMM,
    [np.array([-WIDTH/2, 0]), np.array([WIDTH/2, HEIGHT])],
    [NX, NY],
    cell_type=CellType.quadrilateral,
)
dim = domain.topology.dim 
WRITER.write_mesh(domain)

# -------------- #
# Function Space #
# -------------- #
degree = 2
shape = (dim,)
V = dfx.fem.functionspace(domain, ("P", degree, shape))
u = dfx.fem.Function(V, name="Displacement")
f = dfx.fem.Function(V, name="Force")

# ------------------- #
# Boundary conditions #
# ------------------- #
tol = 1e-14

def clamped_boundary(x):
    return abs(x[1]) < tol


def neumann_boundary(x):
    """
    determines whether a node is on the coupling boundary

    """
    return np.logical_or((np.abs(x[1] - HEIGHT) < tol) , np.abs(np.abs(x[0]) - WIDTH / 2) < tol)


fixed_boundary = dfx.fem.locate_dofs_geometrical(V, clamped_boundary)
coupling_boundary = dfx.mesh.locate_entities_boundary(domain, dim - 1, neumann_boundary)
coupling_boundary_tags = dfx.mesh.meshtags(domain, dim-1, np.sort(coupling_boundary), 1)

bcs = [dfx.fem.dirichletbc(np.zeros((dim,)), fixed_boundary, V)]

# ------------ #
# PRECICE INIT #
# ------------ #
participant = Adapter(MPI_COMM, PARTICIPANT_CONFIG, domain, material_properties=material)
participant.initialize(V, neumann_boundary)
dt = participant.dt


# ------------------------ #
# linear elastic equations #
# ------------------------ #
    
E = dfx.fem.Constant(domain, E)
nu = dfx.fem.Constant(domain, NU)
rho = dfx.fem.Constant(domain, RHO)

lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)

def epsilon(v):
    return ufl.sym(ufl.grad(v))

def sigma(v):
    return lmbda * ufl.tr(epsilon(v)) * ufl.Identity(dim) + 2 * mu * epsilon(v)


# ------------------- #
# Time discretization #
# ------------------- #
# prev time step
u_old = dfx.fem.Function(V)
v_old = dfx.fem.Function(V)
a_old = dfx.fem.Function(V)

# current time step
a_new = dfx.fem.Function(V)
v_new = dfx.fem.Function(V)

beta = dfx.fem.Constant(domain, BETA_)
gamma = dfx.fem.Constant(domain, GAMMA_)
    
dx = ufl.Measure("dx", domain=domain)
ds = ufl.Measure("ds", subdomain_data=coupling_boundary_tags)

a = 1 / beta / dt**2 * (u - u_old - dt * v_old) + a_old * (1 - 1 / 2 / beta)
a_expr = dfx.fem.Expression(a, V.element.interpolation_points())

v = v_old + dt * ((1 - gamma) * a_old + gamma * a)
v_expr = dfx.fem.Expression(v, V.element.interpolation_points())

# ------------------ #
# mass, a stiffness  #
# ------------------ #
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)


def mass(u, u_):
    return rho * ufl.dot(u, u_) * dx

def stiffness(u, u_):
    return ufl.inner(sigma(u), epsilon(u_)) * dx

Residual = mass(a, u_) + stiffness(u, u_) - ufl.dot(f, u_) * ds(1)

Residual_du = ufl.replace(Residual, {u: du})
a_form = ufl.lhs(Residual_du)
L_form = ufl.rhs(Residual_du)

opts = PETSc.Options()  # type: ignore
opts["ksp_type"] = "cg"
opts["ksp_rtol"] = 1.0e-8
opts["pc_type"] = "gamg"

# Use Chebyshev smoothing for multigrid
opts["mg_levels_ksp_type"] = "chebyshev"
opts["mg_levels_pc_type"] = "jacobi"

# Improve estimate of eigenvalues for Chebyshev smoothing
opts["mg_levels_ksp_chebyshev_esteig_steps"] = 10

problem = LinearProblem(a_form, L_form, u=u, bcs=bcs, petsc_options=opts.getAll())


# parameters for Time-Stepping
t = 0.0
n = 0

while participant.is_coupling_ongoing():

    if participant.requires_writing_checkpoint():  # write checkpoint
        participant.store_checkpoint(u_old, v_old, a_old, t)

    read_data = participant.read_data(dt)

    f.vector[participant.interface_dof] = read_data
    f.vector.assemble()
    problem.solve()

    f.x.scatter_forward()
    u.x.scatter_forward()

    # Write new displacements to preCICE
    participant.write_data(u)

    # Call to advance coupling, also returns the optimum time step value
    participant.advance(dt)

    # Either revert to old step if timestep has not converged or move to next timestep
    if participant.requires_reading_checkpoint():
        u_cp, v_cp, a_cp, t_cp, = participant.retrieve_checkpoint()
        u_old.interpolate(u_cp)
        v_old.interpolate(v_cp)
        a_old.interpolate(a_cp)
        t = t_cp
    else:
        v_new.interpolate(v_expr)
        a_new.interpolate(a_expr)
        u.vector.copy(u_old.vector)
        v_new.vector.copy(v_old.vector)
        a_new.vector.copy(a_old.vector)
        t += dt

    if participant.is_time_window_complete():
        if n % 10 == 0:
            WRITER.write_function(u, t)
            WRITER.write_function(f, t)
            WRITER.close()

WRITER.close()
participant.finalize()

That gives the results below. This chart seems to be very good but the problem is that the DOLFINx results are scaled 1:100. But as It follows the same pattern of the deformation, I think that the problem is on the DOLFINx problem implementation, and not in the coupling interface.

image

I posted this question in the preCICE discourse, to have an opinion of the preCICE community about the results. And I received an answer that I can’t use the Linear form in that way, because, the incoming load from the fluid participant is applied locally over DoFs, not over the surface so no need for integration, which leads me to mix both examples, into this one:

import numpy as np
import ufl
import os
from mpi4py import MPI

import dolfinx as dfx
import dolfinx.fem.petsc
from dolfinx.mesh import create_rectangle, CellType
from petsc4py import PETSc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, set_bc
from dolfinx.fem.petsc import LinearProblem
from turbine.solvers.interface import Adapter
from turbine.materials import Material
# --------- #
# CONSTANTS #
# --------- #

MPI_COMM = MPI.COMM_WORLD
CURRENT_FOLDER = os.path.dirname(__file__)
PARTICIPANT_CONFIG = os.path.join(CURRENT_FOLDER, "precice-config.json")
RESULTS_DIR = os.path.join(CURRENT_FOLDER, "results")
os.makedirs(RESULTS_DIR, exist_ok=True)
os.chdir(CURRENT_FOLDER)

WRITER = dfx.io.VTKFile(MPI_COMM, f"{RESULTS_DIR}/result.pvd", "w")

WIDTH, HEIGHT = 0.1, 1 
NX, NY = 4, 26


E = 4000000.0
NU = 0.3
RHO = 3000.0
material =  Material("Flap",E=E,nu=NU,rho=RHO)

BETA_ = 0.25
GAMMA_ = 0.5

# ------- #
# MESHING #
# ------- #

domain = create_rectangle(
    MPI_COMM,
    [np.array([-WIDTH/2, 0]), np.array([WIDTH/2, HEIGHT])],
    [NX, NY],
    cell_type=CellType.quadrilateral,
)
dim = domain.topology.dim 
WRITER.write_mesh(domain)

# -------------- #
# Function Space #
# -------------- #
degree = 2
shape = (dim,)
V = dfx.fem.functionspace(domain, ("P", degree, shape))
u = dfx.fem.Function(V, name="Displacement")
f = dfx.fem.Function(V, name="Force")

# ------------------- #
# Boundary conditions #
# ------------------- #
tol = 1e-14

def clamped_boundary(x):
    return abs(x[1]) < tol


def neumann_boundary(x):
    """
    determines whether a node is on the coupling boundary

    """
    return np.logical_or((np.abs(x[1] - HEIGHT) < tol) , np.abs(np.abs(x[0]) - WIDTH / 2) < tol)


fixed_boundary = dfx.fem.locate_dofs_geometrical(V, clamped_boundary)
coupling_boundary = dfx.mesh.locate_entities_boundary(domain, dim - 1, neumann_boundary)
coupling_boundary_tags = dfx.mesh.meshtags(domain, dim-1, np.sort(coupling_boundary), 1)

bcs = [dfx.fem.dirichletbc(np.zeros((dim,)), fixed_boundary, V)]

# ------------ #
# PRECICE INIT #
# ------------ #
participant = Adapter(MPI_COMM, PARTICIPANT_CONFIG, domain, material_properties=material)
participant.initialize(V, neumann_boundary)
dt = participant.dt


# ------------------------ #
# linear elastic equations #
# ------------------------ #
    
E = dfx.fem.Constant(domain, E)
nu = dfx.fem.Constant(domain, NU)
rho = dfx.fem.Constant(domain, RHO)

lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)

def epsilon(v):
    return ufl.sym(ufl.grad(v))

def sigma(v):
    return lmbda * ufl.tr(epsilon(v)) * ufl.Identity(dim) + 2 * mu * epsilon(v)


# ------------------- #
# Time discretization #
# ------------------- #
# prev time step
u_old = dfx.fem.Function(V)
v_old = dfx.fem.Function(V)
a_old = dfx.fem.Function(V)

# current time step
a_new = dfx.fem.Function(V)
v_new = dfx.fem.Function(V)

beta = dfx.fem.Constant(domain, BETA_)
gamma = dfx.fem.Constant(domain, GAMMA_)
    
dx = ufl.Measure("dx", domain=domain)
ds = ufl.Measure("ds", subdomain_data=coupling_boundary_tags)

a = 1 / beta / dt**2 * (u - u_old - dt * v_old) + a_old * (1 - 1 / 2 / beta)
a_expr = dfx.fem.Expression(a, V.element.interpolation_points())

v = v_old + dt * ((1 - gamma) * a_old + gamma * a)
v_expr = dfx.fem.Expression(v, V.element.interpolation_points())

# ------------------ #
# mass, a stiffness  #
# ------------------ #
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)


def mass(u, u_):
    return rho * ufl.dot(u, u_) * dx

def stiffness(u, u_):
    return ufl.inner(sigma(u), epsilon(u_)) * dx

Residual = mass(a, u_) + stiffness(u, u_) - ufl.dot(f, u_) * ds(1)

Residual_du = ufl.replace(Residual, {u: du})
a_form = ufl.lhs(Residual_du)
a = dfx.fem.form(a_form)

A = assemble_matrix(a, bcs=bcs)
A.assemble()

solver = PETSc.KSP().create(domain.comm)  # type: ignore
solver.setFromOptions()
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.CG)
pc = solver.getPC()
pc.setType(PETSc.PC.Type.GAMG)
solver.setTolerances(rtol=1e-5, atol=1e-11, max_it=300)

# parameters for Time-Stepping
t = 0.0
n = 0

# Resolver el sistema lineal
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setFromOptions()

while participant.is_coupling_ongoing():

    if participant.requires_writing_checkpoint():  # write checkpoint
        participant.store_checkpoint(u_old, v_old, a_old, t)

    read_data = participant.read_data(dt)

    f.x.array[:] = 0
    f.vector[participant.interface_dof] = read_data
    f.vector.assemble()
    set_bc(f.vector, bcs)

    solver.setConvergenceTest(lambda ksp, n, rnorm: print(f"Iteración {n}, Residuo {rnorm}"))
    solver.solve(f.vector, u.vector)
    u.x.scatter_forward()

    # Write new displacements to preCICE
    participant.write_data(u)

    # Call to advance coupling, also returns the optimum time step value
    participant.advance(dt)

    # Either revert to old step if timestep has not converged or move to next timestep
    if participant.requires_reading_checkpoint():
        u_cp, v_cp, a_cp, t_cp, = participant.retrieve_checkpoint()
        u_old.interpolate(u_cp)
        v_old.interpolate(v_cp)
        a_old.interpolate(a_cp)
        t = t_cp
    else:
        v_new.interpolate(v_expr)
        a_new.interpolate(a_expr)
        u.vector.copy(u_old.vector)
        v_new.vector.copy(v_old.vector)
        a_new.vector.copy(a_old.vector)
        t += dt

    if participant.is_time_window_complete():
        if n % 10 == 0:
            WRITER.write_function(u, t)
            WRITER.write_function(f, t)
            WRITER.close()

WRITER.close()
participant.finalize()

But the results of this other version are worse at least the previous one copies the shape of the deformation at another scale, in this case, the solution is completely different in magnitude and shape, and I don’t know where to follow I have tested all that I found here on the discourse forum, and tutorials. I think the main difference between my first code example and my last (from the DOLFINx point of view) is that the first example is a static problem and the second is dynamic, and I don’t know if my mistakes are on how I assemble de linear form in the solve loop.

The only thing that I think that I’m doing is translating the legacy Fenics tutorial to the current version of DOLFINx, using my custom DOLFINx-preCICE adapter, whose main part is the following code that returns the coupling interface DoFs, and cords.

def dof_points_in_vector_space(V, fn):
    bs = V.dofmap.bs

    fs_coords = V.tabulate_dof_coordinates()
    boundary_dofs = fem.locate_dofs_geometrical(V, fn)
    unrolled_dofs = np.empty(len(boundary_dofs) * bs, dtype=np.int32)

    for i, dof in enumerate(boundary_dofs):
        for b in range(bs):
            unrolled_dofs[i * bs + b] = dof * bs + b

    return unrolled_dofs.reshape(-1, bs), fs_coords[boundary_dofs]
1 Like
2

Point-sources can be found in: Point Sources (Redux) - #6 by dokken

3

I solve it by extending the linear solver, like this:

class DiscreteLinearProblem(LinearProblem):
    def __init__(
        self,
        a: ufl.Form,
        L: ufl.Form,
        bcs: list[DirichletBC] = [],
        u: typing.Optional[_Function] = None,
        point_dofs: typing.Optional[dict] = None,
        petsc_options: typing.Optional[dict] = None,
        form_compiler_options: typing.Optional[dict] = None,
        jit_options: typing.Optional[dict] = None,
    ):
        if point_dofs is not None:
            if not isinstance(type(point_dofs), np.ndarray):
                point_dofs = np.array([point_dofs], dtype="int32")
        self._point_dofs = point_dofs
        super().__init__(
            a, L, bcs, u, petsc_options, form_compiler_options, jit_options
        )
    def solve(self, values: typing.Optional[DofsValues] = None):
        """Solve the problem."""
        # Assemble lhs
        self._A.zeroEntries()
        assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
        self._A.assemble()
        # Assemble rhs
        with self._b.localForm() as b_loc:
            b_loc.set(0)
        if self._point_dofs is not None and values is not None:
            # apply load in dofs
            self._b[self._point_dofs] = values
            self._b.assemble()
        assemble_vector(self._b, self._L)
        # Apply boundary conditions to the rhs
        apply_lifting(self._b, [self._a], bcs=[self.bcs])
        self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(self._b, self.bcs)
        # Solve linear system and update ghost values in the solution
        self._solver.solve(self._b, self._x)
        self.u.x.scatter_forward()
        return self.u