Problem with dolfinx_mpc when assemble vetors in paraller

dolfinx
,
1

Hi,

I’m working with the dolfinx and dolfinx_mpc, in CPP interface, from the Poisson demo, I learn that we should assemble the vector as follows:

    b.set(0.0);
    fem::assemble_vector(b.mutable_array(), *L);
    fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, 1.0);
    b.scatter_rev(std::plus<T>());
    fem::set_bc(b.mutable_array(), {bc});

then I want to impose periodic BC, since there is no demo for dolfinx_mpc in CPP, so I tried to write my own code, and I do something like this:

 b.set(0.0);
 dolfinx_mpc::assemble_vector(b.mutable_array(), *L, mpc);
 dolfinx_mpc::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, 1.0, mpc);
 b.scatter_rev(std::plus<T>());
 fem::set_bc(b.mutable_array(), {bc});

The main difference is the second and third lines, the problem is that: when I run my code in 1 process, it works very well, the solution is as expected, but if I try to work in parallel, the error is reported:

[1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[1]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[1]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[1]PETSC ERROR: to get more information on the crash.
[1]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 1 (rank 1 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 1

and this never happen when PBC is not applied. So I’m really confused, could anyone tell me is this part of code corret?
Thanks in advance!

2

Without seeing how you define b there is no way to help you.

3

Hi, thanks for replying, my vector b is defined as follows:

la::Vector<T> b(L->function_spaces()[0]->dofmap()->index_map,
                L->function_spaces()[0]->dofmap()->index_map_bs());

I also find that this error doestn’t happen whe my RHS is 0, that is to say, when all the functions in the form L is 0. So maybe there is some conflict when add the constraint.

4

Hi,

I just did some more test, it seems like it is the problem with the PBC. My computation domain is a 2D rectangle, if I impose PBC in one direction, then everything is good, but when I do it in both dimensions, the error is reported.
I define the PBC as follows:

// PBC_x_condition
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] = (xmax + xmin) - x[i];
    y[i + len] = x[i + len];
    y[i + 2 * len] = x[i + 2 * len];
  }

  return y;
}

and

// PBC_xy_condition
std::vector<std::int8_t> periodic_condition(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) - xmin) < 1e-6) or (std::abs(x(1, p) - ymin) < 1e-6));
  }
  return f;
}

// PBC_xy_relation
std::vector<double> periodic_relation(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++)
  {
    if ((std::abs(x[i] - xmin) < 1e-6) and (std::abs(x[i + len] - ymin) < 1e-6))
    {
      y[i] = (xmax + xmin) - x[i];
      y[i + len] = (ymax + ymin) - x[i + len];
      y[i + 2 * len] = x[i + 2 * len];
    }
    else if (std::abs(x[i] - xmin) < 1e-6)
    {
      y[i] = (xmax + xmin) - x[i];
      y[i + len] = x[i + len];
      y[i + 2 * len] = x[i + 2 * len];
    }
    else if (std::abs(x[i + len] - ymin) < 1e-6)
    {
      y[i] = x[i];
      y[i + len] = (ymax + ymin) - x[i + len];
      y[i + 2 * len] = x[i + 2 * len];
    }
  }

  return y;
}

is there anything wrong with the 2D-PBC? Thanks very much!

5

This is the problem.
DOLFINx_MPC extends the dofmap.
Use the dofmap from the FunctionSpace that is attached to the MPC, ie replace L->function_spaces()[0] with

6

Thanks! Now it works.