Periodic boundary mapping does not work well for some specific settings

dolfinx
1

Hello, I’m working with the dolfinx_mpc and I’m facing some strange phenomenon about the boundary mapping.

Generally I want to creat a 1st-order CG space on a box: [0, 4pi]* [-5,5]* [-5,5], with the grid number Nx, Ny, Nz.
And I want to impose PBC on the x direction, i.e., f(0,y,z) = f(4pi,y,z). So I think the number of slaves node and master node should be equal to (Ny+1)(Nz+1). Then I define my PBC as follows:

const double pi = 3.141592653589793238462;
const double xmin = 0.0;
const double xmax = 4*pi;
const double ymin = -5.0;
const double ymax = 5.0;
const double zmin = -5.0;
const double zmax = 5.0; 
const std::size_t Nx = 10;
const std::size_t Ny = 10;
const std::size_t Nz = 10;

std::array<double,3> p0{xmin,ymin,zmin};
std::array<double,3> p1{xmax,ymax,zmax};

// Create mesh and function space
auto part = mesh::create_cell_partitioner(mesh::GhostMode::none);
auto mesh = std::make_shared<dolfinx::mesh::Mesh>(mesh::create_box(MPI_COMM_WORLD, {{p0, p1}}, {Nx, Ny, Nz}, dolfinx::mesh::CellType::hexahedron, part));
// Define some CG space
V = std::make_shared<fem::FunctionSpace>(fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));

auto pbc = dolfinx_mpc::create_periodic_condition_geometrical(V, periodic_condition_x, periodic_relation_x, {}, 1, false);
auto mpc = std::make_shared<dolfinx_mpc::MultiPointConstraint<T>>(dolfinx_mpc::MultiPointConstraint<double>(V, {pbc.slaves}, {pbc.masters}, {pbc.coeffs}, {pbc.owners}, {pbc.offsets}));
 std::cout<<pbc.slaves.size()<<"\n";

where

std::vector<std::int8_t> periodic_condition_x(std::experimental::mdspan<const double, std::experimental::extents<std::size_t, 3, std::experimental::dynamic_extent>> x)
{
  std::vector<std::int8_t> f;
  for (std::size_t p = 0; p < x.extent(1); ++p)
  {
    f.push_back((std::abs(x(0, p) - xmax) < 1e-8));
  }
  return f;
}

// PBC_x_relation
std::vector<double> periodic_relation_x(std::span<const double> x)
{
  std::vector<double> y(x.size(), 0.0);

  std::size_t len = x.size() / 3;
  for (std::size_t i = 0; i < len; i++)
  {
    y[i] = (xmax + xmin) - x[i];
    y[i + len] = x[i + len];
    y[i + 2 * len] = x[i + 2 * len];
  }

  return y;
}

The problem is that, this code works only for some specifc Nx. In the code above, we know that the number of slaves node is 121, and it should not change when Nx is changed. But if I simply change Nx to 100, then I get 21 slaves dofs. I think I’m not defining the periodic relation in the correct way, does anyone know what could be the problem?
Thanks very much!

2

It would be good to make the problem reproducible by supplying the full C++ code (as you are using the C++ interface).

Another thing you could do for the sake of debugging, is to print what goes into

i.e std::cout<< x[i] << ", " << x[i+len] ", " << x[i+2*len] << " maps to " << y[i] << ", " << y[i+len] ", " << y[i+2*len] << std::endl;
as you would then see what coordinates is found by periodic_condition_x and what they are mapped to in periodic_relation_x.

3

Thanks for reply!
I’m not sure how could I provide the code, should I upload it somewhere? Or I can just paste it here.
I just tested as you said, and it turn out that inside the periodic_boundary_relation, the coordinates of all the correct master dofs and slaves dofs are printed, but after I create the PBC some of them just miss, so I guess some of the slaves dofs could not be found in the function space?

You can play with the parameter “Nx”, by set it to 10 or 100, and you will see that, when Nx = 100, I will not be able to find all the slave-master dofs, while it does when Nx = 10.

Since I didn’t find any demo in C++ interface, this is just my understanding to the source codes, maybe I implemented it in the wrong way.
I have included all of my code here, and please let me know if there is anything more I can supply, thanks!

#include "poisson.h"
#include <PoissonSolver.h>
#include <dolfinx_mpc/PeriodicConstraint.h>

using namespace dolfinx;
using T = PetscScalar;

int main(int argc, char* argv[])
{
  dolfinx::init_logging(argc, argv);
  PetscInitialize(&argc, &argv, nullptr, nullptr);
  {
    std::array<double,3> p0{xmin,ymin,zmin};
    std::array<double,3> p1{xmax,ymax,zmax};

    // Create mesh and function space
    auto part = mesh::create_cell_partitioner(mesh::GhostMode::none);
    auto mesh = std::make_shared<dolfinx::mesh::Mesh>(mesh::create_box(MPI_COMM_WORLD, {{p0, p1}}, {Nx, Ny, Nz}, dolfinx::mesh::CellType::hexahedron, part));

    auto V = std::make_shared<fem::FunctionSpace>(fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));
    
    std::vector<double> coordinates_ = V->tabulate_dof_coordinates(false);
    auto pbc = dolfinx_mpc::create_periodic_condition_geometrical(V, periodic_condition_x, periodic_relation_x, {}, 1, false);
    auto mpc = std::make_shared<dolfinx_mpc::MultiPointConstraint<T>>(dolfinx_mpc::MultiPointConstraint<double>(V, {pbc.slaves}, {pbc.masters}, {pbc.coeffs}, {pbc.owners}, {pbc.offsets}));
   
    std::cout<<pbc.slaves.size()<<"\n";
    for(std::size_t i=0;i<pbc.slaves.size();i++)
    {
std::cout<<coordinates_[3*pbc.slaves[i]]<<" "<<coordinates_[3*pbc.slaves[i]+1]<<" "<<coordinates_[3*pbc.slaves[i]+2]<<" "<<coordinates_[3*pbc.masters[i]]<<" "<<coordinates_[3*pbc.masters[i]+1]<<" "<<coordinates_[3*pbc.masters[i]+2]<<"\n";
    }
   
  }

  PetscFinalize();

  return 0;
}

the source file and python file for forms

#include <dolfinx.h>
using namespace dolfinx;
using T = PetscScalar;


const double pi = 3.141592653589793238462;
const double xmin = 1.0;
const double xmax = 4.0 * pi;
const double ymin = -5.0;
const double ymax = 5.0;
const double zmin = -5.0;
const double zmax = 5.0; 
const std::size_t Nx = 100;
const std::size_t Ny = 2;
const std::size_t Nz = 2;

std::vector<std::int8_t> periodic_condition_x(std::experimental::mdspan<const double, std::experimental::extents<std::size_t, 3, std::experimental::dynamic_extent>> x)
{
  std::vector<std::int8_t> f;
  for (std::size_t p = 0; p < x.extent(1); ++p)
  {
    f.push_back((std::abs(x(0, p) - xmax) < 1e-6));
  }
  return f;
}

// PBC_x_relation
std::vector<double> periodic_relation_x(std::span<const double> x)
{
  std::vector<double> y(x.size(), 0.0);

  std::size_t len = x.size() / 3;
  for (std::size_t i = 0; i < len; i++)
  {
    y[i] = x[i] - (xmax - xmin);
    y[i + len] = x[i + len];
    y[i + 2 * len] = x[i + 2 * len];
    std::cout<<x[i]<<" "<<x[i+len]<<" "<<x[i+2*len]<<" "<<y[i]<<" "<<y[i+len]<<" "<<y[i+2*len]<<std::endl;
  }

  return y;
}

import basix
from basix.ufl_wrapper import create_element
from ufl import Coefficient, TrialFunction, TestFunction, Constant, dx, dS, dot, jump, grad, FacetNormal, dS

cell = "hexahedron"
element = create_element("Lagrange", cell, 1, basix.LagrangeVariant.equispaced)

u  = TrialFunction(element)
w  = TestFunction(element)

f = Coefficient(element)
g = Coefficient(element)
k = Constant(cell)

a = u * w * dx 
L = -g * f.dx(0) * w * dx #- f.dx(1) * w * dx - f.dx(2) * w * dx #- 0.0375 * 0.5 * (f.dx(0) * w.dx(0)+f.dx(1) * w.dx(1)+f.dx(2) * w.dx(2)) * dx 

#------------------------------------------------
4

I just tested the same model with the Python interface (I also include the code), and this problem does not exist, so I think his is not a bug.
Is it possible that this is a issue with the data type?

import dolfinx.fem as fem
import numpy as np
from dolfinx.common import Timer
from dolfinx.mesh import create_box
from dolfinx_mpc import MultiPointConstraint
from mpi4py import MPI
from petsc4py import PETSc

# Get PETSc int and scalar types
complex_mode = True if np.dtype(PETSc.ScalarType).kind == 'c' else False

# Create mesh and finite element
NX = 10
NY = 2
NZ = 2
#mesh = create_unit_cube(MPI.COMM_WORLD, NX, NY, NZ)
mesh = create_box(MPI.COMM_WORLD,((0.0, -5.0, -5.0), (4.0 * np.pi, 5.0, 5.0)),(NX, NY, NZ))
V = fem.FunctionSpace(mesh, ("CG", 1))


bcs = []


def periodic_boundary(x):
    return np.isclose(x[0], 4.0 * np.pi)


def periodic_relation(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0]-4.0 * np.pi
    out_x[1] = x[1]
    out_x[2] = x[2]
    return out_x


with Timer("~PERIODIC: Initialize MPC"):
    mpc = MultiPointConstraint(V)
    mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcs)
    mpc.finalize()


print(mpc.slaves)
5

I have a similar problem here, but the problem does also exist with the Python interface for large meshes, i.e. Nx=Ny=Nz=200.
If you find a solution or the reasion for the problem, please let me know!

1 Like
6

Some update:

It turns out that it is the problem of the round-off error (at least in my case), I can give a little example.

Let’s say I have a cube with the following boundary point

const double xmin = 0.0;
const double xmax = 1.0;
const double ymin = 0.0;
const double ymax = 1.0;
const double zmin = 0.0;
const double zmax = 1.0;

and the PBC relation (x-direction) should be defined in such way

for (std::size_t i = 0; i < len; i++)
  {
    y[i] = (xmax+xmin) - x[i];   # This should be the general definition for PBC
    y[i + len] = x[i + len];
    y[i + 2 * len] = x[i + 2 * len];
  }

I try to map the right bound to left, so in my case
xmax + xmin - x[i] = 1.0 + 0.0 -1.0 = 0.0, but the results turns out to be some tiny number (3.46945e-16) instead of 0.0, that is why MPC can not find the corresponding slaves in the function space.

But the following codes seem to work:

for (std::size_t i = 0; i < len; i++)
  {
    y[i] = xmin;
    y[i + len] = x[i + len];
    y[i + 2 * len] = x[i + 2 * len];
  }

Not sure will it work for you, but I guess it should be the problem about the tolerance or round-off stuff if you have similar problem. Good luck!

2 Likes
7

Thank you for your explanations!
In Python, wouldn’t it be something like:

def pbc_slave_to_master_map_x(x):
    out_x = x.copy()
    out_x[0] = np.round(x[0] - domain.geometry.x[:, 0].max(), decimals=10)
    out_x[1] = x[1]
    out_x[2] = x[2]
    return out_x

Because when I integrate these roundings into my code, the results do not change.

8

There are probably some tolerances that should be adjusted in:

as tol is dolfinx_mpc/cpp/PeriodicConstraint.h at afa9fd51eb7deb8098469fcb3b9185c4ce3099c0 · jorgensd/dolfinx_mpc · GitHub
which might be too small.

1 Like
9

I used higher tolerances until const U tol = 500000 * std::numeric_limits<U>::epsilon(); and it didn’t solve the issue (at least in my case) - the results are the same. Are there more tolerances that might need to be changed?

10

You can search for numerics_limits<U> in the C++ code and see if there are any more that you hit.

11

Does this happen when you use a sparse mesh? like 5 * 5 * 5. For me it is not like a problem of MPC, but the problem with the mesh and function space, in my case, the coordinates of the dofs are not exactly the same as it should be.

12

No, this behavior only occurs for meshes with a large number of elements. For a small mesh (let’s say a unit cube of 50x50x50), the periodic constraints are set correctly. @dokken also suggested in the other post that it is probably related to some rounding issues/too small tolerances, but I couldn’t find the solution for it yet.

13

If your mesh is just cube, maybe you can use out_x[0] = xmin directly, instead of doing this calculation.

14

Unfortunately, this does not work. But thanks for your help!