Error in NewFrac FEniCSx tutorial : ValueError: too many values to unpack (expected 2)

[violetus]

Hello,
I am trying to implement this tutorial related to fracture mechanics and phase field damage models using FEniCSx but while implementing this cell:

load = 0.6
nonzero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])

with nonzero_u.vector.localForm() as bc_local:
    bc_local.set(load) 

problem_u.solve()

state = {"u": u, "alpha": alpha}

plot_damage_state(state, load=load)

I am getting an error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-143-3b01748be6ff> in <module>
      9 state = {"u": u, "alpha": alpha}
     10 
---> 11 plot_damage_state(state, load=load)

~/Documents/FEM_FEniCS_tuts/newfrac-fenicsx-training/notebooks/python/plots.py in plot_damage_state(state, load)
     15     plotter = pyvista.Plotter(title="Damage state", window_size=[800, 300], shape=(1, 2))
     16 
---> 17 # Create VTK unstructured grid from mesh
     18     vtk_mesh = dolfinx.plot.create_vtk_topology(mesh)
     19     grid = pyvista.UnstructuredGrid(vtk_mesh[0], vtk_mesh[1], mesh.geometry.x)

ValueError: too many values to unpack (expected 2)

I tried to fix the file plots.py which calls plot_damage_state function but no luck in suceeding, this is the code in plots.py file by the way:

import numpy as np
import dolfinx.plot
import pyvista

def plot_damage_state(state, load=None):
    """
    Plot the displacement and damage field with pyvista
    """
    u = state["u"]
    alpha = state["alpha"]

    mesh = u.function_space.mesh
    V = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 1))

    plotter = pyvista.Plotter(title="Damage state", window_size=[800, 300], shape=(1, 2))

# Create VTK unstructured grid from mesh
    vtk_mesh = dolfinx.plot.create_vtk_topology(mesh)
    grid = pyvista.UnstructuredGrid(vtk_mesh[0], vtk_mesh[1], mesh.geometry.x)


    plotter.subplot(0, 0)
    if load is not None:
        plotter.add_text(f"Displacement - load {load:3.3f}", font_size=11)
    else:
        plotter.add_text("Displacement", font_size=11)
    vals_2D = u.compute_point_values().real
    vals = np.zeros((vals_2D.shape[0], 3))
    vals[:, :2] = vals_2D
    grid["u"] = vals
    warped = grid.warp_by_vector("u", factor=0.1)
    actor_1 = plotter.add_mesh(warped, show_edges=False)
    plotter.view_xy()

    plotter.subplot(0, 1)
    if load is not None:
        plotter.add_text(f"Damage - load {load:3.3f}", font_size=11)
    else:
        plotter.add_text("Damage", font_size=11)
    grid.point_arrays["alpha"] = alpha.compute_point_values().real
    grid.set_active_scalars("alpha")
    plotter.add_mesh(grid, show_edges=False, show_scalar_bar=True, clim=[0, 1])
    plotter.view_xy()
    pyvista.OFF_SCREEN = False
    plotter.show()

I humbly request to kindly guide me how to solve this issue.
Thank You

[dokken]

Next time, please provide a complete code.
This required me to sit down and convert large portions of the code myself.

import matplotlib.pyplot as plt
import numpy as np

import dolfinx
import dolfinx.plot
import ufl

from mpi4py import MPI
from petsc4py import PETSc

import pyvista
from pyvista.utilities.xvfb import start_xvfb
start_xvfb(wait=0.5)


from plots import plot_damage_state
L = 1.; H = 0.3;
ell_ = 0.1
cell_size = ell_/6;

nx = int(L/cell_size)
ny = int(H/cell_size)

mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (L, H)], [nx, ny])


ndim = mesh.geometry.dim


pyvista.OFF_SCREEN 
dolfinx.plot.create_vtk_mesh(mesh)
grid = pyvista.UnstructuredGrid(*dolfinx.plot.create_vtk_mesh(mesh))
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.add_axes()
plotter.set_scale(5,5)
#plotter.reset_camera(render=True, bounds=(-L/2, L/2, -H/2, H/2, 0, 0))
# if not pyvista.OFF_SCREEN:
#     plotter.show()

# from pathlib import Path
# Path("output").mkdir(parents=True, exist_ok=True)
# figure = plotter.screenshot("output/mesh.png")

element_u = ufl.VectorElement('Lagrange',mesh.ufl_cell(),degree=1,dim=2)
V_u = dolfinx.fem.FunctionSpace(mesh, element_u)

element_alpha = ufl.FiniteElement('Lagrange',mesh.ufl_cell(),degree=1)
V_alpha = dolfinx.fem.FunctionSpace(mesh, element_alpha)

# Define the state
u = dolfinx.fem.Function(V_u, name="Displacement")
alpha = dolfinx.fem.Function(V_alpha, name="Damage")

state = {"u": u, "alpha": alpha}

# need upper/lower bound for the damage field
alpha_lb = dolfinx.fem.Function(V_alpha, name="Lower bound")
alpha_ub = dolfinx.fem.Function(V_alpha, name="Upper bound")

# Measures
dx = ufl.Measure("dx",domain=mesh)
ds = ufl.Measure("ds",domain=mesh)

def bottom(x):
    return np.isclose(x[1], 0.0)

def right(x):
    return np.isclose(x[0], L)

def top(x):
    return np.isclose(x[1], H)

def left(x):
    return np.isclose(x[0], 0.0)

blocked_dofs_left_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), left)
blocked_dofs_right_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), right)
blocked_dofs_bottom_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(1), V_u.sub(1).collapse()[0]), bottom)
blocked_dofs_left_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, left)
blocked_dofs_right_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, right)

zero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])

with zero_u.vector.localForm() as bc_local:
    bc_local.set(0.0)

nonzero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])
with nonzero_u.vector.localForm() as bc_local:
    bc_local.set(1.0)
    
one_alpha = dolfinx.fem.Function(V_alpha)
with one_alpha.vector.localForm() as bc_local:
    bc_local.set(1.0)
                 
zero_alpha = dolfinx.fem.Function(V_alpha)
with zero_alpha.vector.localForm() as bc_local:
    bc_local.set(0.0)


bc_u0 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_left_u, V_u.sub(1))
bc_u1 = dolfinx.fem.dirichletbc(nonzero_u, blocked_dofs_right_u, V_u.sub(1))
bc_u2 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_bottom_u, V_u.sub(0))

bc_alpha0 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_left_alpha)
bc_alpha1 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_right_alpha)

bcs_u = [bc_u0,bc_u1,bc_u2]
bcs_alpha = [bc_alpha0,bc_alpha1]

# setting the upper bound to 0 where BCs are applied
alpha_ub.interpolate(one_alpha)
dolfinx.fem.petsc.set_bc(alpha_ub.vector, bcs_alpha)

E, nu = dolfinx.fem.Constant(mesh, 100.0), dolfinx.fem.Constant(mesh, 0.3)
Gc = dolfinx.fem.Constant(mesh, 1.0)
ell = dolfinx.fem.Constant(mesh, ell_)

def w(alpha):
    """Dissipated energy function as a function of the damage """
    return alpha

def a(alpha, k_ell=1.e-6):
    """Stiffness modulation as a function of the damage """
    return (1 - alpha) ** 2 + k_ell

def eps(u):
    """Strain tensor as a function of the displacement"""
    return ufl.sym(ufl.grad(u))

def sigma_0(u):
    """Stress tensor of the undamaged material as a function of the displacement"""
    mu    = E / (2.0 * (1.0 + nu))
    lmbda = E * nu / (1.0 - nu ** 2)
    return 2.0 * mu * eps(u) + lmbda * ufl.tr(eps(u)) * ufl.Identity(ndim)

def sigma(u,alpha):
    """Stress tensor of the damaged material as a function of the displacement and the damage"""
    return a(alpha) * sigma_0(u)

import sympy 
z = sympy.Symbol("z")
c_w = 4*sympy.integrate(sympy.sqrt(w(z)),(z,0,1))
print("c_w = ",c_w)

c_1w = sympy.integrate(sympy.sqrt(1/w(z)),(z,0,1))
print("c_1/w = ",c_1w)

tmp = 2*(sympy.diff(w(z),z)/sympy.diff(1/a(z),z)).subs({"z":0})
sigma_c = sympy.sqrt(tmp * Gc.value * E.value / (c_w * ell.value))
print("sigma_c = %2.3f"%sigma_c)

eps_c = float(sigma_c/E.value)
print("eps_c = %2.3f"%eps_c)

f = dolfinx.fem.Constant(mesh,(0.,0.))
elastic_energy = 0.5 * ufl.inner(sigma(u,alpha), eps(u)) * dx 
dissipated_energy = Gc / float(c_w) * (w(alpha) / ell + ell * ufl.dot(ufl.grad(alpha), ufl.grad(alpha))) * dx
external_work = ufl.dot(f, u) * dx 
total_energy = elastic_energy + dissipated_energy - external_work   

E_u = ufl.derivative(total_energy,u,ufl.TestFunction(V_u))
E_du = ufl.replace(E_u,{u: ufl.TrialFunction(V_u)})

problem_u = dolfinx.fem.petsc.LinearProblem(a=ufl.lhs(E_du), L=ufl.rhs(E_du), bcs=bcs_u, u=u,
                                      petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

problem_u.solve()

load = .6
with nonzero_u.vector.localForm() as bc_local:
        bc_local.set(load) 
        
problem_u.solve()

plot_damage_state(state,load=load)

and `plots.py

import numpy as np
import dolfinx.plot
import pyvista

def plot_damage_state(state, load=None):
    """
    Plot the displacement and damage field with pyvista
    """
    u = state["u"]
    alpha = state["alpha"]

    mesh = u.function_space.mesh
    pyvista.OFF_SCREEN=True
    plotter = pyvista.Plotter(title="Damage state", window_size=[800, 300], shape=(1, 2))

    # Create VTK unstructured grid from mesh
    vtk_mesh = dolfinx.plot.create_vtk_mesh(u.function_space)
    pyvista.start_xvfb(1)
    grid = pyvista.UnstructuredGrid(*vtk_mesh)


    plotter.subplot(0, 0)
    if load is not None:
        plotter.add_text(f"Displacement - load {load:3.3f}", font_size=11)
    else:
        plotter.add_text("Displacement", font_size=11)
    vals_2D = u.x.array.reshape(-1, u.function_space.dofmap.bs)
    vals = np.zeros((vals_2D.shape[0], 3))
    vals[:, :2] = vals_2D
    grid["u"] = vals
    warped = grid.warp_by_vector("u", factor=0.1)
    actor_1 = plotter.add_mesh(warped, show_edges=False)
    plotter.view_xy()

    plotter.subplot(0, 1)
    if load is not None:
        plotter.add_text(f"Damage - load {load:3.3f}", font_size=11)
    else:
        plotter.add_text("Damage", font_size=11)
    alpha_mesh = pyvista.UnstructuredGrid(*dolfinx.plot.create_vtk_mesh(alpha.function_space))

    alpha_mesh.point_data["alpha"] = alpha.x.array.real
    alpha_mesh.set_active_scalars("alpha")
    plotter.add_mesh(alpha_mesh, show_edges=False, show_scalar_bar=True, clim=[0, 1])
    plotter.view_xy()
    plotter.screenshot("fig.png")

[violetus]

@dokken

Next time, please provide a complete code.

I apologise for not supplying the complete code, next time I will be more careful regarding this.
Thanks for helping me, code is working now but when I tried to implement the next cell which is -

E_alpha = ufl.derivative(total_energy,alpha,ufl.TestFunction(V_alpha))
E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))

damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)

b = dolfinx.la.create_petsc_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
J = dolfinx.fem.form(damage_problem.a)

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

# We set the bound (Note: they are passed as reference and not as values)
solver_alpha_snes.setVariableBounds(alpha_lb.vector,alpha_ub.vector)

I am getting the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-6ac853039d33> in <module>
     11 solver_alpha_snes.setType("vinewtonrsls")
     12 solver_alpha_snes.setFunction(damage_problem.F, b)
---> 13 solver_alpha_snes.setJacobian(damage_problem.J, J)
     14 solver_alpha_snes.setTolerances(rtol=1.0e-9, max_it=50)
     15 solver_alpha_snes.getKSP().setType("preonly")

TypeError: Argument 'J' has incorrect type (expected petsc4py.PETSc.Mat, got Form)

I also tried changing snes_problem.py, but the error still persists, this is my snes_problem.py

from typing import List

from petsc4py import PETSc

import dolfinx
import ufl

        
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, bcs: List[dolfinx.fem.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.form(self.a)
        self.b = dolfinx.fem.form(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.vector)
        self.u.vector.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()

This is the code I am trying to implement:

import matplotlib.pyplot as plt
import numpy as np

import dolfinx
import dolfinx.plot
import ufl

from mpi4py import MPI
from petsc4py import PETSc

import pyvista
from pyvista.utilities.xvfb import start_xvfb
start_xvfb(wait=0.5)

import sys
sys.path.append("/home/ratan/Documents/FEM_FEniCS_tuts/newfrac-fenicsx-training/notebooks/python/")
from snes_problem import SNESProblem

from plots import plot_damage_state
L = 1.; H = 0.3;
ell_ = 0.1
cell_size = ell_/6;

nx = int(L/cell_size)
ny = int(H/cell_size)

mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (L, H)], [nx, ny])


ndim = mesh.geometry.dim


pyvista.OFF_SCREEN 
dolfinx.plot.create_vtk_mesh(mesh, ndim)
grid = pyvista.UnstructuredGrid(*dolfinx.plot.create_vtk_mesh(mesh, ndim))
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.add_axes()
plotter.set_scale(5,5)

#plotter.reset_camera(render=True, bounds=(-L/2, L/2, -H/2, H/2, 0, 0))
# if not pyvista.OFF_SCREEN:
#     plotter.show()

# from pathlib import Path
# Path("output").mkdir(parents=True, exist_ok=True)
# figure = plotter.screenshot("output/mesh.png")

element_u = ufl.VectorElement('Lagrange',mesh.ufl_cell(),degree=1,dim=2)
V_u = dolfinx.fem.FunctionSpace(mesh, element_u)

element_alpha = ufl.FiniteElement('Lagrange',mesh.ufl_cell(),degree=1)
V_alpha = dolfinx.fem.FunctionSpace(mesh, element_alpha)

# Define the state
u = dolfinx.fem.Function(V_u, name="Displacement")
alpha = dolfinx.fem.Function(V_alpha, name="Damage")

state = {"u": u, "alpha": alpha}

# need upper/lower bound for the damage field
alpha_lb = dolfinx.fem.Function(V_alpha, name="Lower bound")
alpha_ub = dolfinx.fem.Function(V_alpha, name="Upper bound")

# Measures
dx = ufl.Measure("dx",domain=mesh)
ds = ufl.Measure("ds",domain=mesh)

def bottom(x):
    return np.isclose(x[1], 0.0)

def right(x):
    return np.isclose(x[0], L)

def top(x):
    return np.isclose(x[1], H)

def left(x):
    return np.isclose(x[0], 0.0)

blocked_dofs_left_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), left)
blocked_dofs_right_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), right)
blocked_dofs_bottom_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(1), V_u.sub(1).collapse()[0]), bottom)
blocked_dofs_left_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, left)
blocked_dofs_right_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, right)

zero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])

with zero_u.vector.localForm() as bc_local:
    bc_local.set(0.0)

nonzero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])
with nonzero_u.vector.localForm() as bc_local:
    bc_local.set(1.0)
    
one_alpha = dolfinx.fem.Function(V_alpha)
with one_alpha.vector.localForm() as bc_local:
    bc_local.set(1.0)
                 
zero_alpha = dolfinx.fem.Function(V_alpha)
with zero_alpha.vector.localForm() as bc_local:
    bc_local.set(0.0)


bc_u0 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_left_u, V_u.sub(1))
bc_u1 = dolfinx.fem.dirichletbc(nonzero_u, blocked_dofs_right_u, V_u.sub(1))
bc_u2 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_bottom_u, V_u.sub(0))

bc_alpha0 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_left_alpha)
bc_alpha1 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_right_alpha)

bcs_u = [bc_u0,bc_u1,bc_u2]
bcs_alpha = [bc_alpha0,bc_alpha1]

# setting the upper bound to 0 where BCs are applied
alpha_ub.interpolate(one_alpha)
dolfinx.fem.petsc.set_bc(alpha_ub.vector, bcs_alpha)

E, nu = dolfinx.fem.Constant(mesh, 100.0), dolfinx.fem.Constant(mesh, 0.3)
Gc = dolfinx.fem.Constant(mesh, 1.0)
ell = dolfinx.fem.Constant(mesh, ell_)

def w(alpha):
    """Dissipated energy function as a function of the damage """
    return alpha

def a(alpha, k_ell=1.e-6):
    """Stiffness modulation as a function of the damage """
    return (1 - alpha) ** 2 + k_ell

def eps(u):
    """Strain tensor as a function of the displacement"""
    return ufl.sym(ufl.grad(u))

def sigma_0(u):
    """Stress tensor of the undamaged material as a function of the displacement"""
    mu    = E / (2.0 * (1.0 + nu))
    lmbda = E * nu / (1.0 - nu ** 2)
    return 2.0 * mu * eps(u) + lmbda * ufl.tr(eps(u)) * ufl.Identity(ndim)

def sigma(u,alpha):
    """Stress tensor of the damaged material as a function of the displacement and the damage"""
    return a(alpha) * sigma_0(u)

import sympy 
z = sympy.Symbol("z")
c_w = 4*sympy.integrate(sympy.sqrt(w(z)),(z,0,1))
print("c_w = ",c_w)

c_1w = sympy.integrate(sympy.sqrt(1/w(z)),(z,0,1))
print("c_1/w = ",c_1w)

tmp = 2*(sympy.diff(w(z),z)/sympy.diff(1/a(z),z)).subs({"z":0})
sigma_c = sympy.sqrt(tmp * Gc.value * E.value / (c_w * ell.value))
print("sigma_c = %2.3f"%sigma_c)

eps_c = float(sigma_c/E.value)
print("eps_c = %2.3f"%eps_c)

f = dolfinx.fem.Constant(mesh,(0.,0.))
elastic_energy = 0.5 * ufl.inner(sigma(u,alpha), eps(u)) * dx 
dissipated_energy = Gc / float(c_w) * (w(alpha) / ell + ell * ufl.dot(ufl.grad(alpha), ufl.grad(alpha))) * dx
external_work = ufl.dot(f, u) * dx 
total_energy = elastic_energy + dissipated_energy - external_work   

E_u = ufl.derivative(total_energy,u,ufl.TestFunction(V_u))
E_du = ufl.replace(E_u,{u: ufl.TrialFunction(V_u)})

problem_u = dolfinx.fem.petsc.LinearProblem(a=ufl.lhs(E_du), L=ufl.rhs(E_du), bcs=bcs_u, u=u,
                                      petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

problem_u.solve()

load = .6
with nonzero_u.vector.localForm() as bc_local:
        bc_local.set(load) 
        
problem_u.solve()

plot_damage_state(state,load=load)

E_alpha = ufl.derivative(total_energy,alpha,ufl.TestFunction(V_alpha))
E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))

damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)

b = dolfinx.la.create_petsc_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
J = dolfinx.fem.form(damage_problem.a)

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

# We set the bound (Note: they are passed as reference and not as values)
solver_alpha_snes.setVariableBounds(alpha_lb.vector,alpha_ub.vector)

I appreciate the guidance to solve this issue.
Thank You for your time.

[dokken]

You have diverged from the original code:
as it states:

J = dolfinx.fem.create_matrix(damage_problem.a)

Ref notebooks/04-phase-field/phase-field.ipynb · main · Newfrac ITN project / newfrac-fenicsx-training · GitLab
while you use

@violetus

J = dolfinx.fem.form(damage_problem.a)

[violetus]

Thank you @dokken for the reference, i followed the link and changed my snes_problem.py file and code accordingly but now this error comes:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-21-e028e09a1a81> in <module>
      2 E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))
      3 
----> 4 damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)
      5 
      6 b = dolfinx.cpp.la.create_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)

~/Documents/FEM_FEniCS_tuts/newfrac-fenicsx-training/notebooks/python/snes_problem.py in __init__(self, F, J, u, bcs)
     28         # Create matrix and vector to be used for assembly
     29         # of the non-linear problem
---> 30         self.A = dolfinx.fem.create_matrix(self.a)
     31         self.b = dolfinx.fem.create_vector(self.L)
     32 

/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in form(form, dtype, form_compiler_params, jit_params)
    137         return form
    138 
--> 139     return _create_form(form)
    140 
    141 

/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in _create_form(form)
    132         return form argument"""
    133         if isinstance(form, ufl.Form):
--> 134             return _form(form)
    135         elif isinstance(form, collections.abc.Iterable):
    136             return list(map(lambda sub_form: _create_form(sub_form), form))

/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in _form(form)
    126                       _cpp.fem.IntegralType.vertex: subdomains.get("vertex")}
    127 
--> 128         return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, code)
    129 
    130     def _create_form(form):

/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py in __init__(self, form, V, coeffs, constants, subdomains, mesh, code)
     46         self._ufcx_form = form
     47         ffi = cffi.FFI()
---> 48         super().__init__(ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
     49                          V, coeffs, constants, subdomains, mesh)
     50 

RuntimeError: Incompatible mesh

Code in snes_problem.py:

from typing import List

from petsc4py import PETSc

import dolfinx
import ufl

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

    def __init__(self, F: ufl.form.Form, J: ufl.form.Form, u: dolfinx.Function, bcs: List[dolfinx.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.create_matrix(self.a)
        self.b = dolfinx.fem.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.vector)
        self.u.vector.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()

Code now -

import matplotlib.pyplot as plt
import numpy as np

import dolfinx
import dolfinx.plot
import ufl

from mpi4py import MPI
from petsc4py import PETSc

import pyvista
from pyvista.utilities.xvfb import start_xvfb
start_xvfb(wait=0.5)

import sys
sys.path.append("/home/ratan/Documents/FEM_FEniCS_tuts/newfrac-fenicsx-training/notebooks/python/")
from snes_problem import SNESProblem

from plots import plot_damage_state
L = 1.; H = 0.3;
ell_ = 0.1
cell_size = ell_/6;

nx = int(L/cell_size)
ny = int(H/cell_size)

mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (L, H)], [nx, ny])


ndim = mesh.geometry.dim


pyvista.OFF_SCREEN 
dolfinx.plot.create_vtk_mesh(mesh, ndim)
grid = pyvista.UnstructuredGrid(*dolfinx.plot.create_vtk_mesh(mesh, ndim))
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.add_axes()
plotter.set_scale(5,5)

#plotter.reset_camera(render=True, bounds=(-L/2, L/2, -H/2, H/2, 0, 0))
# if not pyvista.OFF_SCREEN:
#     plotter.show()

# from pathlib import Path
# Path("output").mkdir(parents=True, exist_ok=True)
# figure = plotter.screenshot("output/mesh.png")

element_u = ufl.VectorElement('Lagrange',mesh.ufl_cell(),degree=1,dim=2)
V_u = dolfinx.fem.FunctionSpace(mesh, element_u)

element_alpha = ufl.FiniteElement('Lagrange',mesh.ufl_cell(),degree=1)
V_alpha = dolfinx.fem.FunctionSpace(mesh, element_alpha)

# Define the state
u = dolfinx.fem.Function(V_u, name="Displacement")
alpha = dolfinx.fem.Function(V_alpha, name="Damage")

state = {"u": u, "alpha": alpha}

# need upper/lower bound for the damage field
alpha_lb = dolfinx.fem.Function(V_alpha, name="Lower bound")
alpha_ub = dolfinx.fem.Function(V_alpha, name="Upper bound")

# Measures
dx = ufl.Measure("dx",domain=mesh)
ds = ufl.Measure("ds",domain=mesh)

def bottom(x):
    return np.isclose(x[1], 0.0)

def right(x):
    return np.isclose(x[0], L)

def top(x):
    return np.isclose(x[1], H)

def left(x):
    return np.isclose(x[0], 0.0)

blocked_dofs_left_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), left)
blocked_dofs_right_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()[0]), right)
blocked_dofs_bottom_u = dolfinx.fem.locate_dofs_geometrical((V_u.sub(1), V_u.sub(1).collapse()[0]), bottom)
blocked_dofs_left_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, left)
blocked_dofs_right_alpha = dolfinx.fem.locate_dofs_geometrical(V_alpha, right)

zero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])

with zero_u.vector.localForm() as bc_local:
    bc_local.set(0.0)

nonzero_u = dolfinx.fem.Function(V_u.sub(0).collapse()[0])
with nonzero_u.vector.localForm() as bc_local:
    bc_local.set(1.0)
    
one_alpha = dolfinx.fem.Function(V_alpha)
with one_alpha.vector.localForm() as bc_local:
    bc_local.set(1.0)
                 
zero_alpha = dolfinx.fem.Function(V_alpha)
with zero_alpha.vector.localForm() as bc_local:
    bc_local.set(0.0)


bc_u0 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_left_u, V_u.sub(1))
bc_u1 = dolfinx.fem.dirichletbc(nonzero_u, blocked_dofs_right_u, V_u.sub(1))
bc_u2 = dolfinx.fem.dirichletbc(zero_u, blocked_dofs_bottom_u, V_u.sub(0))

bc_alpha0 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_left_alpha)
bc_alpha1 = dolfinx.fem.dirichletbc(zero_alpha, blocked_dofs_right_alpha)

bcs_u = [bc_u0,bc_u1,bc_u2]
bcs_alpha = [bc_alpha0,bc_alpha1]

# setting the upper bound to 0 where BCs are applied
alpha_ub.interpolate(one_alpha)
dolfinx.fem.petsc.set_bc(alpha_ub.vector, bcs_alpha)

E, nu = dolfinx.fem.Constant(mesh, 100.0), dolfinx.fem.Constant(mesh, 0.3)
Gc = dolfinx.fem.Constant(mesh, 1.0)
ell = dolfinx.fem.Constant(mesh, ell_)

def w(alpha):
    """Dissipated energy function as a function of the damage """
    return alpha

def a(alpha, k_ell=1.e-6):
    """Stiffness modulation as a function of the damage """
    return (1 - alpha) ** 2 + k_ell

def eps(u):
    """Strain tensor as a function of the displacement"""
    return ufl.sym(ufl.grad(u))

def sigma_0(u):
    """Stress tensor of the undamaged material as a function of the displacement"""
    mu    = E / (2.0 * (1.0 + nu))
    lmbda = E * nu / (1.0 - nu ** 2)
    return 2.0 * mu * eps(u) + lmbda * ufl.tr(eps(u)) * ufl.Identity(ndim)

def sigma(u,alpha):
    """Stress tensor of the damaged material as a function of the displacement and the damage"""
    return a(alpha) * sigma_0(u)

import sympy 
z = sympy.Symbol("z")
c_w = 4*sympy.integrate(sympy.sqrt(w(z)),(z,0,1))
print("c_w = ",c_w)

c_1w = sympy.integrate(sympy.sqrt(1/w(z)),(z,0,1))
print("c_1/w = ",c_1w)

tmp = 2*(sympy.diff(w(z),z)/sympy.diff(1/a(z),z)).subs({"z":0})
sigma_c = sympy.sqrt(tmp * Gc.value * E.value / (c_w * ell.value))
print("sigma_c = %2.3f"%sigma_c)

eps_c = float(sigma_c/E.value)
print("eps_c = %2.3f"%eps_c)

f = dolfinx.fem.Constant(mesh,(0.,0.))
elastic_energy = 0.5 * ufl.inner(sigma(u,alpha), eps(u)) * dx 
dissipated_energy = Gc / float(c_w) * (w(alpha) / ell + ell * ufl.dot(ufl.grad(alpha), ufl.grad(alpha))) * dx
external_work = ufl.dot(f, u) * dx 
total_energy = elastic_energy + dissipated_energy - external_work   

E_u = ufl.derivative(total_energy,u,ufl.TestFunction(V_u))
E_du = ufl.replace(E_u,{u: ufl.TrialFunction(V_u)})

problem_u = dolfinx.fem.petsc.LinearProblem(a=ufl.lhs(E_du), L=ufl.rhs(E_du), bcs=bcs_u, u=u,
                                      petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

problem_u.solve()

load = .6
with nonzero_u.vector.localForm() as bc_local:
        bc_local.set(load) 
        
problem_u.solve()

plot_damage_state(state,load=load)

EE_alpha = ufl.derivative(total_energy,alpha,ufl.TestFunction(V_alpha))
E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))

damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)

b = dolfinx.cpp.la.create_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
J = dolfinx.fem.create_matrix(damage_problem.a)

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

# We set the bound (Note: they are passed as reference and not as values)
solver_alpha_snes.setVariableBounds(alpha_lb.vector,alpha_ub.vector)

Thanks for your patience and time.

[Quentin_Schmid]

Hi,
I tried to implement this example as well. I made some changes to comply with the last version of fenicsx (0.8) but I get an error when trying to solve the problem

here is the code

import matplotlib.pyplot as plt
import numpy as np

import time
import ufl
import sympy 

from dolfinx.mesh import create_rectangle, CellType
from dolfinx.fem import (dirichletbc, Expression, Function, functionspace,\
                         locate_dofs_topological, locate_dofs_geometrical,\
                         set_bc, Constant)
from dolfinx.fem.petsc import LinearProblem    
from dolfinx.io import XDMFFile

from ufl import (TestFunction, TrialFunction, as_vector, dot, dx,\
                 grad, inner, sqrt, split, as_tensor, derivative,\
                 replace, lhs, rhs)

from mpi4py import MPI
from petsc4py import PETSc



import sys
sys.path.append("../python/")
from snes_problem import SNESProblem

t0=time.time()


###############################################################################
########################## MESH ###############################################
###############################################################################
L = 1.; H = 0.3;
ell_ = 0.1
cell_size = ell_/6;

nx = int(L/cell_size)
ny = int(H/cell_size)

msh = create_rectangle(MPI.COMM_WORLD, [(0.0, 0.0), (L, H)], [nx, ny], CellType.triangle)
ndim = msh.geometry.dim

with XDMFFile(MPI.COMM_WORLD, "Mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    
###############################################################################
############################### SPACES ########################################
###############################################################################

V_u = functionspace(msh, ("Lagrange", 1, (msh.geometry.dim, )))
V_alpha = functionspace(msh, ("Lagrange", 1))

u = Function(V_u, name="Displacement")
alpha = Function(V_alpha, name="Damage")

alpha_lb = Function(V_alpha, name="Lower bound")
alpha_ub = Function(V_alpha, name="Upper bound")

# Measures
dx = ufl.Measure("dx",domain=msh)
ds = ufl.Measure("ds",domain=msh)

###############################################################################
########################## BOUNDARY CONDITIONS ################################
###############################################################################

def bottom(x):
    return np.isclose(x[1], 0.0)

def right(x):
    return np.isclose(x[0], L)

def top(x):
    return np.isclose(x[1], H)

def left(x):
    return np.isclose(x[0], 0.0)

blocked_dofs_left_u = locate_dofs_geometrical((V_u.sub(0),V_u), left)
blocked_dofs_right_u = locate_dofs_geometrical((V_u.sub(0),V_u), right)
blocked_dofs_bottom_u = locate_dofs_geometrical((V_u.sub(1),V_u), bottom)
blocked_dofs_left_alpha = locate_dofs_geometrical(V_alpha, left)
blocked_dofs_right_alpha = locate_dofs_geometrical(V_alpha, right)

Vx = V_u.sub(0).collapse()[0]
zero_u = Function(Vx)

with zero_u.vector.localForm() as bc_local:
    bc_local.set(0.0)

nonzero_u = Function(Vx)
with nonzero_u.vector.localForm() as bc_local:
    bc_local.set(1.0)
    
one_alpha = Function(V_alpha)
with one_alpha.vector.localForm() as bc_local:
    bc_local.set(1.0)
                 
zero_alpha = Function(V_alpha)
with zero_alpha.vector.localForm() as bc_local:
    bc_local.set(0.0)


bc_u0 = dirichletbc(zero_u, blocked_dofs_left_u, V_u.sub(1))
bc_u1 = dirichletbc(nonzero_u, blocked_dofs_right_u, V_u.sub(1))
bc_u2 = dirichletbc(zero_u, blocked_dofs_bottom_u, V_u.sub(0))

bc_alpha0 = dirichletbc(zero_alpha, blocked_dofs_left_alpha)
bc_alpha1 = dirichletbc(zero_alpha, blocked_dofs_right_alpha)

bcs_u = [bc_u0,bc_u1,bc_u2]
bcs_alpha = [bc_alpha0,bc_alpha1]

# setting the upper bound to 0 where BCs are applied
alpha_ub.interpolate(one_alpha)
set_bc(alpha_ub.vector, bcs_alpha)

###############################################################################
############################# FORMULATION #####################################
###############################################################################

E, nu = Constant(msh, 100.0), Constant(msh, 0.3)
Gc = Constant(msh, 1.0)
ell = Constant(msh, ell_)

def w(alpha):
    """Dissipated energy function as a function of the damage """
    return alpha

def a(alpha, k_ell=1.e-6):
    """Stiffness modulation as a function of the damage """
    return (1 - alpha) ** 2 + k_ell

def eps(u):
    """Strain tensor as a function of the displacement"""
    return ufl.sym(ufl.grad(u))

def sigma_0(u):
    """Stress tensor of the undamaged material as a function of the displacement"""
    mu    = E / (2.0 * (1.0 + nu))
    lmbda = E * nu / (1.0 - nu ** 2)
    return 2.0 * mu * eps(u) + lmbda * ufl.tr(eps(u)) * ufl.Identity(ndim)

def sigma(u,alpha):
    """Stress tensor of the damaged material as a function of the displacement and the damage"""
    return a(alpha) * sigma_0(u)


z = sympy.Symbol("z")
c_w = 4*sympy.integrate(sympy.sqrt(w(z)),(z,0,1))

c_1w = sympy.integrate(sympy.sqrt(1/w(z)),(z,0,1))

tmp = 2*(sympy.diff(w(z),z)/sympy.diff(1/a(z),z)).subs({"z":0})
sigma_c = sympy.sqrt(tmp * Gc.value * E.value / (c_w * ell.value))

eps_c = float(sigma_c/E.value)

print("c_w = ",c_w)
print("c_1/w = ",c_1w)
print("sigma_c = %2.3f"%sigma_c)
print("eps_c = %2.3f"%eps_c)

f = as_vector([0.0, 0.0])
elastic_energy = 0.5 * inner(sigma(u,alpha), eps(u)) * dx 
dissipated_energy = Gc / float(c_w) * (w(alpha) / ell + ell * dot(grad(alpha), grad(alpha))) * dx
external_work = dot(f, u) * dx 
total_energy = elastic_energy + dissipated_energy - external_work

E_u = derivative(total_energy,u,TestFunction(V_u))
E_du = replace(E_u,{u: TrialFunction(V_u)})

problem_u = LinearProblem(a=lhs(E_du), L=rhs(E_du), bcs=bcs_u, u=u)
problem_u.solve()

#petsc_options={"ksp_type": "preonly", "pc_type": "lu"}

load = .6
with nonzero_u.vector.localForm() as bc_local:
        bc_local.set(load) 
        
problem_u.solve()

with XDMFFile(msh.comm, "U.xdmf", "w") as file:
    file.write_mesh(msh)
    file.write_function(u)

t1=time.time()
t2=t1-t0
h=int(t2/3600)
m=int(t2/60-h*60)
s=t2-m*60-h*3600
print('Overall time = ' + str(h) +' h ' + str(m) + ' min ' + str(s) + ' sec')

which yield the following error

Traceback (most recent call last):
  File "/mnt/c/Users/quent/Documents/Codes/PYTHON_FRESH/FENICSX/PHASE_FIELD/HaleExample.py", line 189, in <module>
    problem_u = LinearProblem(a=lhs(E_du), L=rhs(E_du), bcs=bcs_u, u=u)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/quentin/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 766, in __init__
    self._L = _create_form(
              ^^^^^^^^^^^^^
  File "/home/quentin/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/home/quentin/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/home/quentin/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 178, in _form
    (domain,) = list(sd.keys())  # Assuming single domain
    ^^^^^^^^^
ValueError: not enough values to unpack (expected 1, got 0)
Exception ignored in: <function LinearProblem.__del__ at 0x7ffa9ece54e0>
Traceback (most recent call last):
  File "/home/quentin/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 805, in __del__
    self._solver.destroy()
    ^^^^^^^^^^^^
AttributeError: 'LinearProblem' object has no attribute '_solver'

Any idea on how to solve the problem ?

Thanks,

Quentin