Hello it is me again,
I am trying to something very similar to this
discussion: https://groups.google.com/g/dealii/c/Hh4iBsM21FI/m/cF-mzDP4GkwJ
My question would be how could I pass the solution vector of the
pre-determined solution into the matrix free operator.
I know this is a bad idea but currently I have the vector as a global
variable so I don't have to pass it into the local_apply function or
anything really. I do call initialize_dof_vector with the correct index but
in the conjugate gradient solve function the preconditioner throws the
error included at the bottom of this message.
Also, attached is the just the setup_system() function similar to step-37
to try and keep code to a minimal. Nothing has been altered in the solve()
function from step-37 in the code.
Thanks if anyone has the time to get to this beforehand,
Sean Johnson
--------------------------------------------------------
An error occurred in line <238> of file
</home/sjohnson/AE_beginnings/dealii/include/deal.II/matrix_free/vector_access_internal.h>
in function
void dealii::internal::check_vector_compatibility(const VectorType&,
const dealii::MatrixFree<dim, Number, VectorizedArrayType>&, const
dealii::internal::MatrixFreeFunctions::DoFInfo&) [with int dim = 2; Number
= double; VectorizedArrayType = dealii::VectorizedArray<double, 2>;
VectorType = dealii::LinearAlgebra::distributed::Vector<double>; typename
std::enable_if<has_partitioners_are_compatible<VectorType>,
VectorType>::type* <anonymous> = 0]
The violated condition was:
false
Additional information:
The parallel layout of the given vector is compatible neither with the
Partitioner of the current FEEvaluation with dof_handler_index=1 nor
with any Partitioner in MatrixFree. A potential reason is that you did
not use MatrixFree::initialize_dof_vector() to get a compatible
vector.
Stacktrace:
-----------
#0 ./step-37_minimal: void dealii::internal::check_vector_compatibility<2,
double, dealii::VectorizedArray<double, 2ul>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>,
(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>*)0>(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, dealii::MatrixFree<2, double,
dealii::VectorizedArray<double, 2ul> > const&,
dealii::internal::MatrixFreeFunctions::DoFInfo const&)
#1 ./step-37_minimal: void dealii::FEEvaluationBase<2, 1, double, false,
dealii::VectorizedArray<double, 2ul>
>::read_write_operation<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const, dealii::internal::VectorReader<double,
dealii::VectorizedArray<double, 2ul> >
>(dealii::internal::VectorReader<double, dealii::VectorizedArray<double,
2ul> > const&,
std::array<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const*, 1ul> const&,
std::array<std::vector<dealii::ArrayView<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>::value_type const, dealii::MemorySpace::Host>,
std::allocator<dealii::ArrayView<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>::value_type const, dealii::MemorySpace::Host> >
> const*, 1ul> const&, std::bitset<2ul> const&, bool) const
#2 ./step-37_minimal: void dealii::FEEvaluationBase<2, 1, double, false,
dealii::VectorizedArray<double, 2ul>
>::read_dof_values<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>
>(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, unsigned int, std::bitset<2ul> const&)
#3 ./step-37_minimal: Step37::LaplaceOperator<2, 3,
double>::local_apply(dealii::MatrixFree<2, double,
dealii::VectorizedArray<double, 2ul> > const&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, std::pair<unsigned int, unsigned int>
const&) const
#4 ./step-37_minimal: dealii::internal::MFWorker<dealii::MatrixFree<2,
double, dealii::VectorizedArray<double, 2ul> >,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, Step37::LaplaceOperator<2, 3, double>,
true>::process_range(void (Step37::LaplaceOperator<2, 3, double>::*
const&)(dealii::MatrixFree<2, double, dealii::VectorizedArray<double, 2ul>
> const&, dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, std::pair<unsigned int, unsigned int>
const&) const, std::vector<unsigned int, std::allocator<unsigned int> >
const&, std::vector<unsigned int, std::allocator<unsigned int> > const&,
unsigned int)
#5 ./step-37_minimal: dealii::internal::MFWorker<dealii::MatrixFree<2,
double, dealii::VectorizedArray<double, 2ul> >,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, Step37::LaplaceOperator<2, 3, double>,
true>::cell(unsigned int)
#6 /home/sjohnson/AE_beginnings/dealii/lib/libdeal_II.g.so.9.4.0:
dealii::internal::MatrixFreeFunctions::TaskInfo::loop(dealii::internal::MFWorkerInterface&)
const
#7 ./step-37_minimal: void dealii::MatrixFree<2, double,
dealii::VectorizedArray<double, 2ul>
>::cell_loop<Step37::LaplaceOperator<2, 3, double>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >(void (Step37::LaplaceOperator<2, 3,
double>::*)(dealii::MatrixFree<2, double, dealii::VectorizedArray<double,
2ul> > const&, dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, std::pair<unsigned int, unsigned int>
const&) const, Step37::LaplaceOperator<2, 3, double> const*,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, bool) const
#8 ./step-37_minimal: Step37::LaplaceOperator<2, 3,
double>::apply_add(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&) const
#9 ./step-37_minimal: dealii::MatrixFreeOperators::Base<2,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, dealii::VectorizedArray<double, 2ul>
>::mult_add(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, bool) const
#10 ./step-37_minimal: dealii::MatrixFreeOperators::Base<2,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, dealii::VectorizedArray<double, 2ul>
>::vmult_add(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&) const
#11 ./step-37_minimal: dealii::MatrixFreeOperators::Base<2,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, dealii::VectorizedArray<double, 2ul>
>::vmult(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&) const
#12 ./step-37_minimal:
dealii::internal::SolverCG::IterationWorker<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, Step37::LaplaceOperator<2, 3, double>,
dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >, int>::do_iteration(unsigned int)
#13 ./step-37_minimal: void
dealii::SolverCG<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >::solve<Step37::LaplaceOperator<2, 3, double>,
dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> > >(Step37::LaplaceOperator<2, 3, double>
const&, dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&,
dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> > const&)
#14 ./step-37_minimal:
dealii::PreconditionChebyshev<Step37::LaplaceOperator<2, 3, double>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>,
dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >
>::estimate_eigenvalues(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&) const
#15 ./step-37_minimal:
dealii::PreconditionChebyshev<Step37::LaplaceOperator<2, 3, double>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>,
dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >
>::vmult(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&) const
#16 ./step-37_minimal:
dealii::mg::SmootherRelaxation<dealii::PreconditionChebyshev<Step37::LaplaceOperator<2,
3, double>, dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>,
dealii::DiagonalMatrix<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> > >,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >::apply(unsigned int,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&) const
#17 /home/sjohnson/AE_beginnings/dealii/lib/libdeal_II.g.so.9.4.0:
dealii::Multigrid<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >::level_v_step(unsigned int)
#18 /home/sjohnson/AE_beginnings/dealii/lib/libdeal_II.g.so.9.4.0:
dealii::Multigrid<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >::level_v_step(unsigned int)
#19 /home/sjohnson/AE_beginnings/dealii/lib/libdeal_II.g.so.9.4.0:
dealii::Multigrid<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >::cycle()
#20 ./step-37_minimal: void
dealii::internal::PreconditionMGImplementation::vmult<2,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, dealii::MGTransferMatrixFree<2, double>,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >(std::vector<dealii::DoFHandler<2, 2> const*,
std::allocator<dealii::DoFHandler<2, 2> const*> > const&,
dealii::Multigrid<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> >&, dealii::MGTransferMatrixFree<2, double>
const&, dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&, bool, dealii::mg::Signals const&, ...)
#21 ./step-37_minimal: void dealii::PreconditionMG<2,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, dealii::MGTransferMatrixFree<2, double>
>::vmult<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>
>(dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>&,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host> const&) const
#22 ./step-37_minimal:
dealii::internal::SolverCG::IterationWorker<dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, Step37::LaplaceOperator<2, 3, double>,
dealii::PreconditionMG<2,
dealii::LinearAlgebra::distributed::Vector<double,
dealii::MemorySpace::Host>, dealii::MGTransferMatrixFree<2, double> >,
int>::do_iteration(unsigned int)
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/82328209-1ed1-4760-bcaa-2b818dd395bcn%40googlegroups.com.
template <int dim>
void LaplaceProblem<dim>::setup_system()
{
Timer time;
setup_time = 0;
system_matrix.clear();
mg_matrices.clear_elements();
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
dof_handler_DG.distribute_dofs(fe_DG);
dof_handler_DG.distribute_mg_dofs();
pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
const IndexSet locally_relevant_dofs =
DoFTools::extract_locally_relevant_dofs(dof_handler);
const IndexSet locally_relevant_dofs_DG =
DoFTools::extract_locally_relevant_dofs(dof_handler_DG);
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(
mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
constraints.close();
constraints_DG.clear();
constraints_DG.reinit(locally_relevant_dofs_DG);
DoFTools::make_hanging_node_constraints(dof_handler_DG, constraints_DG);
VectorTools::interpolate_boundary_values(
mapping, dof_handler_DG, 0, Functions::ZeroFunction<dim>(),
constraints_DG);
constraints_DG.close();
setup_time += time.wall_time();
time_details << "Distribute DoFs & B.C. (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << 's' << std::endl;
time.restart();
const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler,
&dof_handler_DG};
const std::vector<const AffineConstraints<double> *> constraints_list =
{&constraints, &constraints_DG};
{
typename MatrixFree<dim, double>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree<dim, double>::AdditionalData::none;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values | update_quadrature_points);
std::shared_ptr<MatrixFree<dim, double>> system_mf_storage(
new MatrixFree<dim, double>());
system_mf_storage->reinit(mapping,
dof_handlers,
constraints_list,
QGauss<1>(fe.degree + 1),
additional_data);
system_mf_storage->initialize_dof_vector(test_solution, 1);
system_matrix.initialize(system_mf_storage,{0});
}
system_matrix.initialize_dof_vector(extra_test);
system_matrix.initialize_dof_vector(solution);
system_matrix.initialize_dof_vector(system_rhs);
system_matrix.evaluate_coefficient(Coefficient<dim>());
setup_time += time.wall_time();
time_details << "Setup matrix-free system (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << 's' << std::endl;
time.restart();
// Next, initialize the matrices for the multigrid method on all the
// levels. The data structure MGConstrainedDoFs keeps information about
// the indices subject to boundary conditions as well as the indices on
// edges between different refinement levels as described in the step-16
// tutorial program. We then go through the levels of the mesh and
// construct the constraints and matrices on each level. These follow
// closely the construction of the system matrix on the original mesh,
// except the slight difference in naming when accessing information on
// the levels rather than the active cells.
const unsigned int nlevels = triangulation.n_global_levels();
mg_matrices.resize(0, nlevels - 1);
const std::set<types::boundary_id> dirichlet_boundary_ids = {0};
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
dirichlet_boundary_ids);
mg_constrained_dofs_DG.initialize(dof_handler_DG);
mg_constrained_dofs_DG.make_zero_boundary_constraints(dof_handler_DG,
dirichlet_boundary_ids);
const std::vector<MGConstrainedDoFs> mg_constrained_list =
{mg_constrained_dofs, mg_constrained_dofs_DG};
for (unsigned int level = 0; level < nlevels; ++level)
{
const IndexSet relevant_dofs =
DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
AffineConstraints<double> level_constraints;
level_constraints.reinit(relevant_dofs);
level_constraints.add_lines(
mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
const IndexSet relevant_dofs_DG =
DoFTools::extract_locally_relevant_level_dofs(dof_handler_DG, level);
AffineConstraints<double> level_constraints_DG;
level_constraints_DG.reinit(relevant_dofs_DG);
level_constraints_DG.add_lines(
mg_constrained_dofs_DG.get_boundary_indices(level));
level_constraints_DG.close();
const std::vector<const AffineConstraints<double> *>
level_constraints_list = {&level_constraints, &level_constraints_DG};
typename MatrixFree<dim, float>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
MatrixFree<dim, float>::AdditionalData::none;
additional_data.mapping_update_flags =
(update_gradients | update_JxW_values | update_quadrature_points);
additional_data.mg_level = level;
std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(
new MatrixFree<dim, float>());
mg_mf_storage_level->reinit(mapping,
dof_handlers,
level_constraints_list,
QGauss<1>(fe.degree + 1),
additional_data);
mg_mf_storage_level->initialize_dof_vector(test_solution,1);
mg_matrices[level].initialize(mg_mf_storage_level,
mg_constrained_dofs,
level,{0});
mg_matrices[level].evaluate_coefficient(Coefficient<dim>());
}
setup_time += time.wall_time();
time_details << "Setup matrix-free levels (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << 's' << std::endl;
}
