DOLFINX_MPC produce different results when work in parallel

Sorry, it’s me agin. It is a continuous question based on my last post, I have now fix the issuse with assembling the vector, now there is not error reported when I run the program.

But I meet another problem: the solver gives me different solution when work in parallel.
Since this doesn’t happend when the MPC is not applied, so I guess there could be some bugs within my codes, maybe there are some ghost nodes need to be updated but I forgot to do that?

Here is the code for creating the matrix and vector

auto pattern = dolfinx_mpc::create_sparsity_pattern(a, mpc, mpc);
auto A = std::make_shared<Matrix>(MPI_COMM_WORLD, pattern);

auto b = std::make_shared<Vector<T>>(mpc->function_space()->dofmap()->index_map,

Then I assemle the matrix and vector, and solve the system as follows:


dolfinx_mpc::assemble_matrix(Matrix::set_block_fn(A.mat(), ADD_VALUES), Matrix::set_fn(A.mat(), ADD_VALUES), *a, mpc, mpc, {}, 1.0);
MatAssemblyBegin(A.mat(), MAT_FLUSH_ASSEMBLY);
MatAssemblyEnd(A.mat(), MAT_FLUSH_ASSEMBLY);

fem::set_diagonal<T>(Matrix::set_fn(A.mat(), INSERT_VALUES), *V, {});
MatAssemblyBegin(A.mat(), MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A.mat(), MAT_FINAL_ASSEMBLY);



dolfinx_mpc::assemble_vector(b.mutable_array(), *L, mpc);
dolfinx_mpc::apply_lifting(b.mutable_array(), {a}, {{}}, {}, 1.0, mpc);

fem::set_bc(b.mutable_array(), {});

la::petsc::Vector _b(la::petsc::create_vector_wrap(*b), false);
la::petsc::Vector _u(la::petsc::create_vector_wrap(*u.x()), false);

lu->solve(_u.vec(), _b.vec());

I would greatly appreciate any comments or suggestions. Thank you very much!

There should be a scatter forward after set_bc.

Also there should be one after solve, Pre-backsubstitute.

Please note that with only snippets, its hard to diagnose issues/bugs in your code

Thanks for replying!

May I ask when should the function scatter_fwd() be applied? I read the Poisson demo here, and it seems to me that, this scatter_fwd function is not applied after set_bc, I guess it is because we only need to scatter forward when work in parallel?

And this doesn’t really solve my issue, the solution is still not correct after paralization.

Thank you very much!

Scatter forward should be applied whenever there has been done any work on only the owned degrees of freedom of a vector (not ghost values), and these have to be updated from the owning processors.

Use-cases where this might be needed:

  1. Interpolation on subset of cells

  2. After scatter_reverse to update accumulated values from owned to ghosted indices.

  3. After using set_bc, as this only updates the owned degrees of freedom.

  4. Using petsc to solve a LA problem (as petsc doesn’t update ghosts after solve)

Without a reproducible example, it is hard to give any further guidance.

Thanks for your patience and suggestions!

Now I include all of the codes, if we consider the simplest case:
f = g.
g = cos(2pi x) in [0, 1] with periodic BC.
Then the solution of f should be a cosine function as well, i.e., we should get f(0) = f(1) = 1

But when I solve this equation in parallel, it seems like the slave node is not really updated, I actually get f(1) = 5.68175e-322, which is close to 0.
Could you please have a look at the codes? Thanks a lot!

using namespace dolfinx;
using T = PetscScalar;

int main(int argc, char *argv[])
  dolfinx::init_logging(argc, argv);
  PetscInitialize(&argc, &argv, nullptr, nullptr);
    // Create mesh and function space
    auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
    auto mesh = std::make_shared<mesh::Mesh>(mesh::create_interval(MPI_COMM_WORLD, 10 , {0, 1.0}));

    auto V = std::make_shared<fem::FunctionSpace>(
        fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));

    auto f = std::make_shared<fem::Function<T>>(V);

    // Define variational forms
    std::map<std::string, std::shared_ptr<const fem::Function<T>>> empty_map;
    auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(*form_poisson_a, {V, V}, empty_map, {}, {}));
    auto L = std::make_shared<fem::Form<T>>(fem::create_form<T>(*form_poisson_L, {V}, {{"f", f}}, {}, {}));

    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}));

    // Define boundary condition

        [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
          std::vector<T> f;
          for (std::size_t p = 0; p < x.extent(1); ++p)
            double dx = x(0, p);

          return {f, {f.size()}};

    // Compute solution
    fem::Function<T> u(V);
    auto pattern = dolfinx_mpc::create_sparsity_pattern(*a,mpc,mpc);

    std::shared_ptr<la::petsc::Matrix> A;
    A = std::make_shared<la::petsc::Matrix>(MPI_COMM_WORLD,pattern);

    la::Vector<T> b(mpc->function_space()->dofmap()->index_map,

    dolfinx_mpc::assemble_matrix(la::petsc::Matrix::set_block_fn(A->mat(), ADD_VALUES), la::petsc::Matrix::set_fn(A->mat(), ADD_VALUES), *a, mpc,mpc, {}, 1.0);
    MatAssemblyBegin(A->mat(), MAT_FLUSH_ASSEMBLY);
    MatAssemblyEnd(A->mat(), MAT_FLUSH_ASSEMBLY);
    fem::set_diagonal<T>(la::petsc::Matrix::set_fn(A->mat(), INSERT_VALUES), *V,
    MatAssemblyBegin(A->mat(), MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A->mat(), MAT_FINAL_ASSEMBLY);

    dolfinx_mpc::assemble_vector(b.mutable_array(), *L, mpc);
    dolfinx_mpc::apply_lifting(b.mutable_array(), {a}, {{}}, {}, 1.0, mpc);
    fem::set_bc(b.mutable_array(), {});

    la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
    la::petsc::options::set("ksp_type", "preonly");
    la::petsc::options::set("pc_type", "lu");

    la::petsc::Vector _u(la::petsc::create_vector_wrap(*u.x()), false);
    la::petsc::Vector _b(la::petsc::create_vector_wrap(b), false);
    lu.solve(_u.vec(), _b.vec());

    VecView(_u.vec(), PETSC_VIEWER_STDOUT_WORLD); 

    for(int i=0;i<u.x()->map()->size_local();i++)
    //  std::cout<<u.x()->array()[i]<<"\n";



  return 0;

and the PBC is defined like this

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;

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

  return y;

Is it possible that your problem is also related to this post? In my case, the slave DOFs are also not updated by the backsubstitution process.

Hi, thanks for sharing! I will read through it.

So have you solved it now?Or do you know how could fix it?
It seems like the slaves dof is not in the same process with the master dof, so it is not properly updated.

As I just posted in the other post, you need to use a function with the MPC function space to store your solution, if you want back-substitution to work properly, for the same reason of using the MPC function space for the RHS vector.


Yes I just found that!
By using u(mpc->function_space()), now it works! Thanks very much!