Hello, thanks for replying! I have read the page that you suggested, and I was able to get some solution of singular system with PBC, but the soltion seems to be very wird,
Is this a bug? or maytbe I have implement something wrong.
I attached my codes here:
int main(int argc, char *argv[])
{
dolfinx::init_logging(argc, argv);
PetscInitialize(&argc, &argv, nullptr, nullptr);
{
// Define basic parameters
auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
auto mesh = std::make_shared<mesh::Mesh>(mesh::create_interval(MPI_COMM_WORLD, N, {xmin, xmax}, part));
// Define function space
auto V = std::make_shared<fem::FunctionSpace>(fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));
// Prepare and set Constants for the bilinear form
auto f = std::make_shared<fem::Function<T>>(V);
auto u = 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}}, {}, {}));
f->interpolate(source);
io::VTKFile file1(MPI_COMM_WORLD, "results/f.pvd", "w");
file1.write<T>({*f}, 0.0);
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}));
LinearProblem Poisson_eq(a, L, mpc, {});
Poisson_eq.assemble_Mat();
Poisson_eq.assemble_Vec();
Poisson_eq.solve(*u);
io::VTKFile file2(MPI_COMM_WORLD, "results/u.pvd", "w");
file2.write<T>({*u}, 0.0);
}
PetscFinalize();
return 0;
}
Here are my source codes:
using namespace dolfinx;
using T = PetscScalar;
//-----------------------------------------------------------------------------
LinearProblem::LinearProblem(
std::shared_ptr<fem::Form<T>> a,
std::shared_ptr<fem::Form<T>> L,
std::shared_ptr<dolfinx_mpc::MultiPointConstraint<T>> mpc,
std::vector<std::shared_ptr<const fem::DirichletBC<T>>> bcs) : _a(a), _L(L), _mpc(mpc), _bcs(bcs),
_b(L->function_spaces()[0]->dofmap()->index_map,
L->function_spaces()[0]->dofmap()->index_map_bs())
{
auto pattern = dolfinx_mpc::create_sparsity_pattern(*a,mpc,mpc);
pattern.assemble();
_matA = std::make_shared<la::petsc::Matrix>(MPI_COMM_WORLD,pattern);
lu = std::make_shared<la::petsc::KrylovSolver>(MPI_COMM_WORLD);
la::petsc::options::set("ksp_type", "gmres");
la::petsc::options::set("pc_type", "hypre");
lu->set_from_options();
}
//-----------------------------------------------------------------------------
LinearProblem::~LinearProblem()
{
}
//-----------------------------------------------------------------------------
void LinearProblem::assemble_Mat()
{
MatZeroEntries(_matA->mat());
dolfinx_mpc::assemble_matrix(la::petsc::Matrix::set_block_fn(_matA->mat(), ADD_VALUES), la::petsc::Matrix::set_fn(_matA->mat(), ADD_VALUES), *_a, _mpc, _mpc, {}, 1.0);
MatAssemblyBegin(_matA->mat(), MAT_FLUSH_ASSEMBLY);
MatAssemblyEnd(_matA->mat(), MAT_FLUSH_ASSEMBLY);
fem::set_diagonal<T>(la::petsc::Matrix::set_fn(_matA->mat(), INSERT_VALUES), *_L->function_spaces()[0], {_bcs});
MatAssemblyBegin(_matA->mat(), MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(_matA->mat(), MAT_FINAL_ASSEMBLY);
lu->set_operator(_matA->mat());
}
//-----------------------------------------------------------------------------
void LinearProblem::assemble_Vec()
{
_b.set(0.0);
dolfinx_mpc::assemble_vector(_b.mutable_array(), *_L, _mpc);
dolfinx_mpc::apply_lifting(_b.mutable_array(), {_a}, {{_bcs}}, {}, 1.0, _mpc);
_b.scatter_rev(std::plus<T>());
fem::set_bc(_b.mutable_array(), {_bcs});
}
//-----------------------------------------------------------------------------
void LinearProblem::solve(fem::Function<T> &u)
{
la::petsc::Vector _bv(la::petsc::create_vector_wrap(_b), false);
la::petsc::Vector _u(la::petsc::create_vector_wrap(*u.x()), false);
PetscBool has_cnst = PETSC_TRUE;
PetscInt num_vectors = 0;
Vec * vec=NULL;
MatNullSpace nullspace;
PetscErrorCode ierr = MatNullSpaceCreate(MPI_COMM_WORLD, has_cnst, num_vectors, vec, &nullspace);
MatSetNullSpace(_matA->mat(),nullspace);
MatNullSpaceRemove(nullspace, _bv.vec());
if (ierr == 0)
MatNullSpaceDestroy(&nullspace);
lu->solve(_u.vec(), _bv.vec());
_mpc->backsubstitution(u.x()->mutable_array());
}
//-----------------------------------------------------------------------------
std::pair<std::vector<T>, std::vector<std::size_t>> source(std::experimental::mdspan<const double, std::experimental::extents<std::size_t, 3, std::experimental::dynamic_extent>> x)
{
std::vector<T> u0;
for (std::size_t p = 0; p < x.extent(1); ++p)
{
double dx = x(0, p);
u0.push_back(4.0 * pi * pi * sin(2.0 * pi * dx));
}
return {u0, {u0.size()}};
}
//-----------------------------------------------------------------------------
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;
}
import basix
from basix.ufl_wrapper import create_element
from ufl import (Coefficient, TestFunction, TrialFunction, dx, ds, inner, grad, interval)
cell = "interval"
degree = 1
element = create_element("Lagrange",cell, degree)
u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
a = inner(grad(u),grad(v)) * dx
L = f*v*dx
Thanks very much!
