Too Many values to unpack (expected 2)

jhssyb · November 7, 2022, 6:34pm

Hi, I tried to get the base flow of the flow past cylinder problems (Re=40). When solving steady-state Naiver-Stokes equations using Taylor-Hood element, I got an error complaining ’ ValueError: too many values to unpack (expected 2) '. I have no idea what is wrong with my code. Hope to get some help from you.

I use Dolfinx v0.5.0. The boundary conditions and mesh generations are exactly the same in the following tutorial as I just copied a snippet from that tutorial.

https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code2.html

However, when I changed the trial function as ‘(u, p) = w.split()’ instead of ‘u, p = TrialFunctions(W)’. The ‘too many value to unpack’ error disappeared. But the solver diverges as all the residuals turn to be NAN. Quite weird, for me, these two statements are almost the same. I could not figure out what is wrong with my solver.

import gmsh
import numpy as np
from dolfinx import nls, log
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from dolfinx import fem
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (Constant, assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,NonlinearProblem,
                               create_vector, create_matrix, set_bc)
from dolfinx.graph import create_adjacencylist
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (distribute_entity_data, gmshio, XDMFFile)
from dolfinx.mesh import create_mesh, meshtags_from_entities

from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunctions, TrialFunctions, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, Identity)


nev = 50 # number of eigen values
Re = 40
nu = 1 / Re


gmsh.initialize()
print('Gmesh version ', gmsh.__version__)
L = 2.2
H = 0.41
c_x = c_y =0.2
r = 0.05
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0

if mesh_comm.rank == model_rank:
    rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1)
    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

if mesh_comm.rank == model_rank:
    fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()
fluid_marker = 1
if mesh_comm.rank == model_rank:
    volumes = gmsh.model.getEntities(dim=gdim)
    assert(len(volumes) == 1)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
    gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    print('boundary ',boundaries)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H/2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H/2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L/2, H, 0]) or np.allclose(center_of_mass, [L/2, 0, 0]):
            walls.append(boundary[1])
        else:
            obstacle.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

# Create distance field from obstacle.
res_min = r / 5
if mesh_comm.rank == model_rank:
    distance_field = gmsh.model.mesh.field.add("Distance")
    gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
    threshold_field = gmsh.model.mesh.field.add("Threshold")
    gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
    gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
    gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H)
    gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
    gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
    min_field = gmsh.model.mesh.field.add("Min")
    gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
    gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

if mesh_comm.rank == model_rank:
    gmsh.option.setNumber("Mesh.Algorithm", 6) # Algorithm 8 will cause MPI_ABORT--PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
    gmsh.option.setNumber("Mesh.RecombineAll", 1)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")

mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
fdim = mesh.topology.dim - 1

# Taylor-Hood element [P2, P2]-P1
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2*P1
# Create the function space
W = fem.FunctionSpace(mesh, TH)
V = fem.FunctionSpace(mesh, P2)
Q = fem.FunctionSpace(mesh, P1)

w = fem.Function(W)
#(u, p) = w.split()
u, p = TrialFunctions(W)
v, q = TestFunctions(W)


# Define boundary conditions
class InletVelocity():
    def __init__(self, velocity):
        self.max_u = velocity

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]),dtype=PETSc.ScalarType)
        values[0] = 4 * self.max_u * x[1] * (0.41 - x[1])/(0.41**2)
        return values

# Inlet
max_velocity = 0.3
u_inlet = fem.Function(V)
inlet_velocity = InletVelocity(max_velocity)
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)),V)

# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V)
bcu = [bcu_inflow, bcu_obstacle, bcu_walls]
# Outlet
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)
bcp = [bcp_outlet]

bcs = [bcu_inflow, bcu_obstacle, bcu_walls]

# Weak form of Steady Navier-Stokes
F = (Constant(mesh,nu) * inner(grad(u), grad(v)) + dot(dot(grad(u), u), v) - p * div(v) + q * div(u)) * dx

# Creating the solver object
problem = NonlinearProblem(F, w, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

# We can customize the linear solver used inside the NewtonSolver by modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(w)
assert(converged)
print(f"Number of interations: {n:d}")

nate · November 7, 2022, 6:42pm

Components of your code don’t make sense. You’re trying to solve it as a nonlinear problem, but you’re attempting to form the bilinear formulation as if it were linear instead of the residual

This causes more problems here when you define your NonlinearProblem

jhssyb:

# Creating the solver object
problem = NonlinearProblem(F, w, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

If you want to solve your system using Newton’s method, consider @dokken’s tutorials: A nonlinear Poisson equation — FEniCSx tutorial

jhssyb · November 7, 2022, 6:48pm

Thanks, nate. This is nonlinear problem and cannot be solved by linear solvers. I formulated this problem as F==0 rather than using the bilinear formula, i.e., a==L. It is recommended to use taylor-hood element to solve steady-state Navier-Stokes, which is different from the tutorial you referred. I am quite confused about the different between these two ways to define ‘u, p’
(u, p) = w.split()
u, p = TrialFunctions(W)

jhssyb · November 7, 2022, 6:51pm

I do not understand why this causes more problems as I just copied this code from the tutorial for nonlinear problems. From my understanding, this is the standard way to formulate and setup nonlinear solvers. Could you point out my misunderstanding?

# Creating the solver object
problem = NonlinearProblem(F, w, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

nate · November 7, 2022, 7:05pm

You edited your post after I commented. Your previous version employed

u, p = TrialFunctions(W)

followed by an attempt to explicitly define a bilinear formulation which naturally cannot be bilinear.

Your corrections allow the Newton solver to start; however, clearly you don’t achieve convergence. It’s also not obvious if your residual is well defined since you immediately generate invalid data in the Newton solver:

[main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = inf (tol = 1e-10) r (rel) = -nan(tol = 1e-06)

There are many thing to check: BCs, formulation, initial guess, well-posed-ness in conjunction with the incompressibility constraint, etc…

dokken · November 7, 2022, 7:09pm

The tutorial uses the Taylor-Hood element, as the velocity space is a second order Lagrange space, and the pressure space is a first Order Lagrange space, see the first lines under: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial

However, is uses a fractional step method, as it is often favoured, as the mixed/coupled non-linear problem is usually quite memory intensive.

jhssyb · November 8, 2022, 8:57am

Thanks, dokken. For unsteady problems, we can use IPCS and calculate velocity and pressure separately using linear solvers in three steps. However, regarding the steady nonlinear problem, it seems to me that we need to solve NS for w and then get the solution using u,v=w.split() ? Can we use the same strategy as IPCS to solve this steady problem?

jhssyb · November 8, 2022, 9:01am

Thanks, nate. How can we setup residual as it is automatically computed from F==0, should we customize residuals? As to u, p = TrialFunctions(W), does this statement perform same function as u,p=w.split(). For me, they are just two different ways to get trial functions, not expecting the first one causing ‘too many values to unpack’ and the other works without such errors.

jhssyb · November 8, 2022, 9:46am

By the way, could you help me check the BCs, I suspect there is sth wrong with this setting but could not identify errors. For this steady problem, no initial guess is required.

dokken · November 8, 2022, 10:01am

I would suggest you have a look at: dolfinx/demo_stokes.py at main · FEniCS/dolfinx · GitHub
which shows you how to set up the variational form for a coupled problem.
You can combine this with:
Implementation — FEniCSx tutorial
or
Custom Newton solvers — FEniCSx tutorial
to create your fully coupled solver.

When you solve non-linear problems, you cannot use TrialFunctions in your variational form. You have to use a function, and use u, p = ufl.split(w), as the residual F has to be differentiated when forming the Jacobian used in Newton’s method.

jhssyb · November 8, 2022, 11:06am

Thanks, I followed the ‘non-linear Poission’ tutorial to create a nonlinear solver. Since this is nonlinear solver using ‘F==0’, I do not think it needs to assemble LHS matrix and RHS vector and then apply lifting. To create a variational form, I use

# Taylor-Hood element [P2, P2]
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2*P1
# Create the function space
W = fem.FunctionSpace(mesh, TH)
V = fem.FunctionSpace(mesh, P2)
Q = fem.FunctionSpace(mesh, P1)

w = fem.Function(W)
u, p = ufl.split(w)
v, q = TestFunctions(W)
# Define variational forms
I = Identity(gdim)
T = 2.0 * Constant(mesh,nu) * sym(grad(u)) - p * I
F = inner(T, grad(v)) * dx - q * div(u) * dx + inner(grad(u) * u, v) * dx

As for the Newton solver, I wrote the following code

# Creating the solver object
problem = NonlinearProblem(F, w, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

# We can customize the linear solver used inside the NewtonSolver by modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(w)
assert(converged)
print(f"Number of interations: {n:d}")

But it does not converge, still producing NAN residual, totally confused

dokken · November 8, 2022, 11:36am

Solving F==0 involves solving J dw = -F(w^k)
Where w^(k+1) = w^k + dw
as I write out explicitly in

The boundary conditions you have supplied above does not make sense, as you Need to consider the BCs in sub spaces, see

jhssyb · November 8, 2022, 11:49am

Thanks, I do not understand why we need to set boundary conditions in subspace. In the unsteady tutorial, boundary conditions in space V,Q were set as

u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)),V)

In order to set it in subspace, you mean we should do this?

bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological((W.sub(0), V), fdim, ft.find(wall_marker)),V)

By the way, why should we customize Newton Solvers, can we use the exact solver as the ‘nonlinear Poisson problem’ tutorial? I am completely new to Dolfinx, sorry for this silly quesiton

dokken · November 8, 2022, 11:52am

As explained in the Stokes demo this should be:

wall_dofs =  locate_dofs_topological((W.sub(0), V), fdim, ft.find(wall_marker))
bcu_walls = dirichletbc(u_nonslip, wall_dofs, W.sub(0))

I am simply trying to give you a notion about what happens under the hood of the non-linear solver, as you stated:

which isn’t correct, as this is done under the hood (and is illustrated in the tutorial I am referencing).

jhssyb · November 8, 2022, 4:42pm

Thanks, a small question about the Stoke tutorial. As to the no slip boundary condition, why do we define noslip in the vector space V using ‘noslip = Function(V)’ but in the second part, the velocity was defined in W0=W.sub(0).collapse() rather than V, as ‘lid_velocity = Function(W0)’. What are the differences between these two statements

# No slip boundary condition
noslip = Function(V)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = dirichletbc(noslip, dofs, W.sub(0))


# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = dirichletbc(lid_velocity, dofs, W.sub(0))

dokken · November 8, 2022, 4:53pm

Those two statements are equivalent in this context.
It is probably a slip-up from the author (should be W0 for consistency) since the demo contains multiple versions of the same problem.

jhssyb · November 8, 2022, 5:08pm

Thanks, I see. I found that you have provided a test file for nonliner assembler. In the class definition of NonlinearPDE_SNESProblem(), what are the purposes of using apply_lifting and ghost_update. After reading the dolfinx document, I still have confusion.

github.com

FEniCS/dolfinx/blob/main/python/test/unit/fem/test_nonlinear_assembler.py

# Copyright (C) 2019 Nathan Sime
#
# This file is part of DOLFINx (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later
"""Unit tests for assembly"""

import math

import numpy as np
import pytest

import ufl
from dolfinx.cpp.la.petsc import scatter_local_vectors
from dolfinx.fem import (Function, FunctionSpace, VectorFunctionSpace,
                         bcs_by_block, dirichletbc, extract_function_spaces,
                         form, locate_dofs_topological)
from dolfinx.fem.petsc import (apply_lifting, apply_lifting_nest,
                               assemble_matrix, assemble_matrix_block,
                               assemble_matrix_nest, assemble_vector,

This file has been truncated. show original

        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)

dokken · November 8, 2022, 5:34pm

You can find the documentation of apply_lifting at dolfinx.fem — DOLFINx 0.5.1 documentation

It is a way of enforcing Dirichlet BCs in such a Wat that your matrix on the LHS is kept symmetric.
As all these operations are done on cells local to the process they are on (when running in parallel with MPI, one has to sum up contributions for dofs that are shared between processes. This is what Ghost update does.
You can read about how dolfinx works with MPI at @jackhale and C. Maurini’s tutorial DOLFINx in Parallel with MPI — NewFrac FEniCSx Training

dg0609 · February 5, 2023, 1:08am

with

W0 = W.sub(0).collapse()

the following line of code

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(W0)

gives the error

AttributeError: 'tuple' object has no attribute '_cpp_object'

That error is avoided with

W0,_ = W.sub(0).collapse()

So is the code above the correct approach?

dokken · February 5, 2023, 7:14am

Collapse returns a tuple. First; the new smaller FunctionSpace. Secondly, a map from each dof in the collapsed space to those in the parent space.

Topic		Replies	Views
ValueError: too many values to unpack (expected 1), while solving 2D elasticity problem Errors	4	511	September 27, 2022
Eigenvalue elasticity problem: "too many values to unpack" error when defining geometric stiffness dolfinx	1	178	December 29, 2023
Mixed Finite element ValueError too many values to unpack dolfinx	2	261	October 17, 2022
Too many values to unpack! General	2	363	December 9, 2022
ValueError: too many values to unpack (expected 2) in solver part Linear Algebra	3	1073	January 7, 2022

Too Many values to unpack (expected 2)

Related topics