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
f->interpolate(
[](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);
f.push_back(cos(2*std::numbers::pi*dx));
}
return {f, {f.size()}};
});
// Compute solution
fem::Function<T> u(V);
auto pattern = dolfinx_mpc::create_sparsity_pattern(*a,mpc,mpc);
pattern.assemble();
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,
mpc->function_space()->dofmap()->index_map_bs());
MatZeroEntries(A->mat());
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);
b.set(0.0);
dolfinx_mpc::assemble_vector(b.mutable_array(), *L, mpc);
dolfinx_mpc::apply_lifting(b.mutable_array(), {a}, {{}}, {}, 1.0, mpc);
b.scatter_rev(std::plus<T>());
fem::set_bc(b.mutable_array(), {});
b.scatter_fwd();
la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
la::petsc::options::set("ksp_type", "preonly");
la::petsc::options::set("pc_type", "lu");
lu.set_from_options();
lu.set_operator(A->mat());
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());
u.x()->scatter_fwd();
mpc->backsubstitution(u.x()->mutable_array());
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";
}
}
PetscFinalize();
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;
}