Petsc error with calling dof_to_vertex

liuyun · July 15, 2021, 4:22pm

I am trying to load a set of data using numpy and converting it to fenics mesh before solving the poisson equation. However, the following error is encountered

Traceback (most recent call last):
File “sample_poisson.py”, line 92, in
coeff0.vector().set_local(coeff0_val[d2v])
RuntimeError:

*** Error: Unable to successfully call PETSc function ‘VecSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /home/yl398/aur/fenics/dolfin/src/dolfin-2019.1.0.post0/dolfin/la/PETScVector.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset: acf365fd0a6e5ea80243a27621f7d9678f01d99f

I have attached a code here. The m_h data contains 13^3=2197 entries, which matches the dimension of the mesh generated for the cube. There should be no problems with the boundary conditions. Any help is greatly appreciated.

from fenics import *
import numpy as np
from math import *

deg = 1
a = 12
## Create mesh and define function space
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(a, a, a), a, a, a)

m_h_input = np.genfromtxt('m_h.dat')


class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool ((near(x[0], 0) or near(x[1], 0) or near(x[2], 0)) and 
            (not ((near(x[0], a) and near(x[2], a)) or 
                  (near(x[0], a) and near(x[1], a)) or
                  (near(x[1], a) and near(x[2], a)))) and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
    	#### define mapping for a single point in the box, such that 3 mappings are required
        if near(x[0], a) and near(x[1], a) and near(x[2],a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2] - a
        ##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
        if near(x[0], a) and near(x[2], a):
            y[0] = x[0] - a
            y[1] = x[1] 
            y[2] = x[2] - a      
        elif near(x[1], a) and near(x[2], a):
            y[0] = x[0] 
            y[1] = x[1] - a
            y[2] = x[2] - a
        elif near(x[0], a) and near(x[1], a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2]         
        #### right maps to left: left/right is defined as the x-direction
        elif near(x[0], a):
            y[0] = x[0] - a
            y[1] = x[1]
            y[2] = x[2]
        ### back maps to front: front/back is defined as the y-direction    
        elif near(x[1], a):
            y[0] = x[0]
            y[1] = x[1] - a
            y[2] = x[2] 
        #### top maps to bottom: top/bottom is defined as the z-direction        
        elif near(x[2], a):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - a

## Define boundary condition
pbc = PeriodicBoundary()

V = FunctionSpace(mesh, 'Lagrange', deg, constrained_domain=pbc)
u = TrialFunction(V)
v = TestFunction(V)

d2v = dof_to_vertex_map(V)
coeff0_val = m_h_input
coeff0 = Function(V)
coeff0.vector().set_local(coeff0_val[d2v])

a = coeff0*inner(grad(u), grad(v))*dx
f = Constant(-6.0)
L = f*v*dx

## Compute solution
u = Function (V)
solve(a == L, u, [])

dokken · July 17, 2021, 11:13am

As you are constraining your function space with a periodic boundary condition, you reduce the number of degrees of freedom in your space. Following your example, you can see this by printing the dim of a constrained and unconstrained function space:

from fenics import *
import numpy as np
from math import *

deg = 1
a = 12
## Create mesh and define function space
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(a, a, a), a, a, a)

m_h_input = np.arange(2197, dtype=float)
#np.genfromtxt('m_h.dat')


class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool ((near(x[0], 0) or near(x[1], 0) or near(x[2], 0)) and 
            (not ((near(x[0], a) and near(x[2], a)) or 
                  (near(x[0], a) and near(x[1], a)) or
                  (near(x[1], a) and near(x[2], a)))) and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
    	#### define mapping for a single point in the box, such that 3 mappings are required
        if near(x[0], a) and near(x[1], a) and near(x[2],a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2] - a
        ##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
        if near(x[0], a) and near(x[2], a):
            y[0] = x[0] - a
            y[1] = x[1] 
            y[2] = x[2] - a      
        elif near(x[1], a) and near(x[2], a):
            y[0] = x[0] 
            y[1] = x[1] - a
            y[2] = x[2] - a
        elif near(x[0], a) and near(x[1], a):
            y[0] = x[0] - a
            y[1] = x[1] - a
            y[2] = x[2]         
        #### right maps to left: left/right is defined as the x-direction
        elif near(x[0], a):
            y[0] = x[0] - a
            y[1] = x[1]
            y[2] = x[2]
        ### back maps to front: front/back is defined as the y-direction    
        elif near(x[1], a):
            y[0] = x[0]
            y[1] = x[1] - a
            y[2] = x[2] 
        #### top maps to bottom: top/bottom is defined as the z-direction        
        elif near(x[2], a):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - a

## Define boundary condition
pbc = PeriodicBoundary()

V_unc = FunctionSpace(mesh, 'Lagrange', deg)
V = FunctionSpace(mesh, 'Lagrange', deg, constrained_domain=pbc)
print(V_unc.dim(), V.dim())

which returns

2197 1797

liuyun · July 19, 2021, 10:44am

Thank you very much for the explanation. That definitely helps a lot. I am having difficulties determining the exact degree of freedom using the periodic boundary conditions.

Initially there were a^3 degrees of freedom, but I’m mapping the 3 faces of the cube to another 3 faces, reducing the degree of freedom by 3a^2. I’m mapping the edges separately which sort of avoids the double counting, so adding an addition 6a degrees of freedom (since there are 12 edges).

So that should be a^3-3a^2+6a = 1768. This is different from the 1797 returned by the code. What am I doing wrong? Is the periodic boundary condition defined correctly?

dokken · July 19, 2021, 12:58pm

I would suggest you start with a slightly simpler mapping, (for instance map one facet to another facet), and try to understand what is going on from that mapping. I do not have time to double check that all your mappings are correct.

However, I think there might be a bug in DOLFIN when it comes to periodic mappings, as the number of dofs are not as expected in this MWE:

from fenics import *
deg = 1
a = 2
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(a, a, a), a, a, a)

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return near(x[0], 0, 1e-6)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1.0
        y[1] = x[1]
        y[2] = x[2]

## Define boundary condition
pbc = PeriodicBoundary()

V_unc = FunctionSpace(mesh, 'Lagrange', deg)
V = FunctionSpace(mesh, 'Lagrange', deg, constrained_domain=pbc)


print(V_unc.dim(), V.dim())
print(V.tabulate_dof_coordinates())
print(V.dim(),(a+1)**3-(a+1)**2, (a+1)**3, V_unc.dim())
assert(V.dim() == (a+1)**3-(a+1)**2)

One thing you can always do to backward engineer the location of the input data is to tabulate dof coordinates (works on CG-1 spaces), and find the matching vertex.

You could also try out my implementation of periodic boundary conditions in dolfinx, namely:

github.com

jorgensd/dolfinx_mpc/blob/master/python/demos/demo_periodic_geometrical.py

# This demo program solves Poisson's equation
#
#     - div grad u(x, y) = f(x, y)
#
# on the unit square with homogeneous Dirichlet boundary conditions
# at y = 0, 1 and periodic boundary conditions at x = 0, 1.
#
# Original implementation in DOLFIN by Kristian B. Oelgaard and Anders Logg
# This implementation can be found at:
# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/periodic/demo_periodic.py
#
# Copyright (C) Jørgen S. Dokken 2020.
#
# This file is part of DOLFINX_MPCX.
#
# SPDX-License-Identifier:    LGPL-3.0-or-later


import dolfinx
import dolfinx.io

This file has been truncated. show original

liuyun · July 26, 2021, 12:58pm

Thank you very much for the suggestions. I will look at the example carefully.

It seems that I would need to install dolfinx for to use the code snippets? Is dolfinx stable for daily use?

dokken · July 26, 2021, 1:02pm

The first example is using dolfin, not dolfinx.

My implementation in dolfinx of periodic conditions are regulary updated to use the latest syntax of dolfinx. We have one release of dolfinx,

Which is matched by:

Topic		Replies	Views
PETSc function error Errors	5	468	April 29, 2022
PETSs Error Code 63 when trying to do Nonlinear Problem	2	684	February 15, 2022
Error: Unable to successfully call PETSc function 'MatSetValuesLocal'. Reason: PETSc error code is: 63 (Argument out of range) Errors	2	2117	July 6, 2021
PETSc error code is: 63 (Argument out of range)	1	1455	April 11, 2019
PETSc error code 63 (argument out of range)	2	3055	December 5, 2019

Petsc error with calling dof_to_vertex

Related Topics