Problems SNESVI converting to newer dolfinx

dolfinx
1

Hi,
I have a phase-field formulation which can reach arbitrary negative and positive values, yet I want it to be constrained to [0, 1]. I looked at different approaches to enforce this, and since pyadjoint isn’t implemented for dolfinx, I believe that the SNES VI solver might be best.
I’m using this tutorial on VI in dolfinx as reference for the solver, but it is giving me several issues in getting it to work with dolfinx 0.8.0. Here are two examples:

  • dolfinx.fem.petsc.create_matrix
AttributeError: 'Form' object has no attribute '_cpp_object'

I’ve seen this post mentioning the need for a fem.form, but I believe that is not the issue.

  • dolfinx.cpp.la.create_vector
AttributeError: module 'dolfinx.cpp.la' has no attribute 'create_vector'

And I expect that there will be more issues if these are resolved.

Here is a MWE, with a commented working unconstrained example, and the constrained example failing.
Any help is much appreciated!

import numpy as np
import dolfinx.fem.petsc
from dolfinx import mesh, fem
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from dolfinx.io import XDMFFile
import ufl
from petsc4py import PETSc
from typing import List

print(f"dolfinx version: {dolfinx.__version__}")

dt = 2e-1
epsilon = 1e-2

mesh_domain = mesh.create_unit_square(MPI.COMM_WORLD, 50, 50, mesh.CellType.triangle)

"""
Function spaces
"""
V = fem.functionspace(mesh_domain, ("CG", 1))         # works as expected

phi = fem.Function(V, name="phase")
phi_prev = fem.Function(V, name="phase_old")
v = ufl.TestFunction(V)

"""
Phase field formulation
"""
phi_transformed = (phi + 0.2)/1.4       # Rescaled phi for the double-well term to make minima at -0.2, 1.2 -> to test whether the constrained optimization keeps phi in [0, 1]
F = (phi - phi_prev) / dt * v * ufl.dx + epsilon**2 * ufl.inner(ufl.grad(phi), ufl.grad(v)) * ufl.dx + phi_transformed * (1 - phi_transformed)*(1 - 2*phi_transformed) * v * ufl.dx

"""
Initial conditions
"""
rng = np.random.default_rng(42)
phi.interpolate(lambda x: 0.5 + 0.5 * (0.5 - rng.random(x.shape[1])))

"""
Unconstrained solver setup:     Working
"""
# problem = NonlinearProblem(F, phi)
# solver = NewtonSolver(MPI.COMM_WORLD, problem)

"""
Constrained solver setup:   Not working
"""
dF = ufl.derivative(F, phi, ufl.TestFunction(V))
ddF = ufl.derivative(dF, phi, ufl.TrialFunction(V))

class SNESProblem:
    """Nonlinear problem class compatible with PETSC.SNES solver."""

    def __init__(self, F: ufl.form.Form, J: ufl.form.Form, u: dolfinx.fem.function.Function, bcs: List[dolfinx.fem.bcs.DirichletBC] = []):
        """This class set up structures for solving a non-linear problem using Newton's method.

        Parameters
        ==========
        F: Residual.
        J: Jacobian.
        u: Solution.
        bcs: Dirichlet boundary conditions.
        """
        self.L = F
        self.a = J
        self.bcs = bcs
        self.u = u

        # Create matrix and vector to be used for assembly
        # of the non-linear problem

        self.A = dolfinx.fem.petsc.create_matrix(self.a)
        self.b = dolfinx.fem.petsc.create_vector(self.L)

    def F(self, snes: PETSc.SNES, x: PETSc.Vec, b: PETSc.Vec):
        """Assemble the residual F into the vector b.

        Parameters
        ==========
        snes: the snes object
        x: Vector containing the latest solution.
        b: Vector to assemble the residual into.
        """
        # We need to assign the vector to the function
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.x)
        self.u.x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        # Zero the residual vector
        with b.localForm() as b_local:
            b_local.set(0.0)
        dolfinx.fem.assemble_vector(b, self.L)

        # Apply boundary conditions
        dolfinx.fem.apply_lifting(b, [self.a], [self.bcs], [x], -1.0)
        b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        dolfinx.fem.set_bc(b, self.bcs, x, -1.0)

    def J(self, snes, x: PETSc.Vec, A: PETSc.Mat, P: PETSc.Mat):
        """Assemble the Jacobian matrix.

        Parameters
        ==========
        x: Vector containing the latest solution.
        A: Matrix to assemble the Jacobian into.
        """
        A.zeroEntries()
        dolfinx.fem.assemble_matrix(A, self.a, self.bcs)
        A.assemble()

snes_problem = SNESProblem(dF, ddF, phi)

b = dolfinx.cpp.la.create_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
J = dolfinx.fem.petsc.create_matrix(snes_problem.a)

# Create Newton solver and solve
solver_snes = PETSc.SNES().create()
solver_snes.setType("vinewtonrsls")
solver_snes.setFunction(snes_problem.F, b)
solver_snes.setJacobian(snes_problem.J, J)
solver_snes.setTolerances(rtol=1.0e-9, max_it=50)
solver_snes.getKSP().setType("preonly")
solver_snes.getKSP().setTolerances(rtol=1.0e-9)
solver_snes.getKSP().getPC().setType("lu")

# We set the bound (Note: they are passed as reference and not as values)
zero = dolfinx.fem.function.Function(V)
with zero.x.localForm() as loc:
    loc.set(0.0)

one = dolfinx.fem.function.Function(V)
with one.x.localForm() as loc:
    loc.set(1.0)
solver_snes.setVariableBounds(zero.x,one.x)

"""
Solving
"""
file = XDMFFile(MPI.COMM_WORLD, "phasefield_constrained.xdmf", "w")
file.write_mesh(mesh_domain)

t = 0.0
timesteps = 50

file.write_function(phi, 0)

for i in range(timesteps):
    phi_prev.x.array[:] = phi.x.array
    t += dt

    # Unconstrained solver
    # r = solver.solve(phi)
    # print(f"Step {int(t / dt)}: num iterations: {r[0]}")

    # Constrained solver
    r = solver_snes.solve(None, phi.x)

    # Process & save solution
    file.write_function(phi, i+1)

file.close()
2

Starting point: you definitely want to use create matrix/vector from dolfinx.fem.petsc, because the one without petsc is creating scipy matrices and numpy arrays.

1 Like
3

Thank you. I tried this already but it results in the same error. I’ve updated the code above to reflect it.

4

You need to wrap the ufl form into a dolfinx.fem.form. See for instance the classes at

1 Like
5

Thank you! Those links were very helpful, I’ve managed to run the problem using parts of it, this is my new code:

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

from dolfinx.fem import form
from dolfinx.fem.petsc import (
    apply_lifting,
    assemble_matrix,
    assemble_vector,
    create_matrix,
    set_bc,
)
from dolfinx.la import create_petsc_vector
from ufl import TrialFunction, derivative

"""
Constrained solver setup:
"""

class NonlinearPDE_SNESProblem:
    def __init__(self, F, u, bc):

        V = u.function_space
        du = TrialFunction(V)
        self.L = form(F)
        self.a = form(derivative(F, u, du))
        self.bc = bc
        self._F, self._J = None, None
        self.u = u

    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.x.petsc_vec)
        self.u.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        assemble_vector(F, self.L)
        apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(F, [self.bc], x, -1.0)

    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        assemble_matrix(J, self.a, bcs=[self.bc])
        J.assemble()

print(f"dolfinx version: {dolfinx.__version__}")

dt = 2e-1
epsilon = 1e-2

mesh_domain = mesh.create_unit_square(MPI.COMM_WORLD, 50, 50, mesh.CellType.triangle)
V = fem.functionspace(mesh_domain, ("CG", 1))
phi = fem.Function(V, name="phase")
phi_prev = fem.Function(V, name="phase_old")
v = ufl.TestFunction(V)

phi_transformed = (phi + 0.2)/1.4       # Rescaled phi for the double-well term to make minima at -0.2, 1.2 -> to test whether the constrained optimization keeps phi in [0, 1]
F = (phi - phi_prev) / dt * v * ufl.dx + epsilon**2 * ufl.inner(ufl.grad(phi), ufl.grad(v)) * ufl.dx + phi_transformed * (1 - phi_transformed)*(1 - 2*phi_transformed) * v * ufl.dx

rng = np.random.default_rng(42)
phi.interpolate(lambda x: 0.5 + 0.5 * (0.5 - rng.random(x.shape[1])))


# Fake BC for testing
u_bc = fem.Function(V)
u_bc.x.array[:] = 1.0
bc = dolfinx.fem.dirichletbc(
    u_bc,
    dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], -1.0)),
)

# Create nonlinear problem
problem = NonlinearPDE_SNESProblem(F, phi, bc)

b = create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
J = create_matrix(problem.a)

# Create Newton solver and solve
solver_snes = PETSc.SNES().create()
solver_snes.setType("vinewtonrsls")
solver_snes.setFunction(problem.F, b)
solver_snes.setJacobian(problem.J, J)
solver_snes.setTolerances(rtol=1.0e-9, max_it=10)
solver_snes.getKSP().setType("preonly")
solver_snes.getKSP().setTolerances(rtol=1.0e-9)
solver_snes.getKSP().getPC().setType("lu")

"""
Solving
"""
file = XDMFFile(MPI.COMM_WORLD, "phasefield_constrained.xdmf", "w")
file.write_mesh(mesh_domain)

t = 0.0
timesteps = 50

file.write_function(phi, 0)

for i in range(timesteps):
    phi_prev.x.array[:] = phi.x.array
    t += dt

    # Constrained solver
    r = solver_snes.solve(None, phi.x.petsc_vec)

    # Process & save solution
    file.write_function(phi, i+1)

file.close()

I have 2 things I haven’t managed to figure out here:

  1. How can I not specify any bcs? In other examples I see that passing an empty list works, however, that leads to an error here:
AttributeError: 'list' object has no attribute '_cpp_object'
  1. I’m having problems with giving the correct type to setVariableBounds:
zero = fem.Function(V)
zero.x.array[:] = 0.0
one = fem.Function(V)
one.x.array[:] = 1.0

solver_snes.setVariableBounds(zero.x,one.x)

Gives:

TypeError: Argument 'xl' has incorrect type (expected petsc4py.PETSc.Vec, got Vector)

I see other examples using something like:

with zero.vector.localForm() as loc:
    loc.set(0.0)

, which gives the same error and uses the deprecated vector.
Alternatively, when I use:

solver_snes.setVariableBounds(petsc4py.PETSc.Vec(zero.x),petsc4py.PETSc.Vec(one.x))

I get a segmentation violation.

6

rewrite this to

class NonlinearPDE_SNESProblem:
    def __init__(self, F, u, bcs):

        V = u.function_space
        du = TrialFunction(V)
        self.L = form(F)
        self.a = form(derivative(F, u, du))
        self.bcs = bcs
        self._F, self._J = None, None
        self.u = u

    def F(self, snes, x, F):
        """Assemble residual vector."""
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        x.copy(self.u.x.petsc_vec)
        self.u.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

        with F.localForm() as f_local:
            f_local.set(0.0)
        assemble_vector(F, self.L)
        apply_lifting(F, [self.a], bcs=[self.bcs], x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(F, self.bcs, x, -1.0)

    def J(self, snes, x, J, P):
        """Assemble Jacobian matrix."""
        J.zeroEntries()
        assemble_matrix(J, self.a, bcs=self.bcs)
        J.assemble()

Now you can pass in an empty list.

Secondly,

This should be zero.x.petsc_vec and one.x.petsc_vec.

1 Like
7

Perfect, that works. Thank you