I have limited the problem down to matrix-free sizes. The reinit() function
seems to tell the matrix-free object to expect to calculate the number of
elements that would fit the larger of the two DoFHandler's handed to the
reinit() function. However, I need to specify that it should actually match
the smaller of the two. Is there a way to accomplish this? I need to pass
both DoFHandlers to the reinit() function because the elements of the
matrix are calculated using variables that are solved from a hyperbolic
equation (thus they have DG elements and a DG DoFHandler) and this specific
equation being solved is an elliptic equation and therefore CG elements are
used and a CG DoFHandler is used.
Please let me know if this is possible or any work arounds someone else has
used.
Thanks,
Sean Johnson
On Thursday, June 22, 2023 at 2:22:27 PM UTC-6 Sean Johnson wrote:
> This is the setup_system() function that has been changed just for the
> system_matrix in step-37.
> dof_handler and fe are the names of the CG version and the DG versions end
> with "_DG"
> system_matrix.clear();
> mg_matrices.clear_elements();
>
> dof_handler.distribute_dofs(fe);
> dof_handler.distribute_mg_dofs();
>
> 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);
> constraints_DG.close();
>
> const std::vector<const DoFHandler<dim> *> dof_handlers =
> {&dof_handler_DG, &dof_handler};
> const std::vector<const AffineConstraints<double> *>
> constraints_list = {&constraints_DG, &constraints};
> 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();
> {
> typename MatrixFree<dim, double>::AdditionalData additional_data;
> additional_data.tasks_parallel_scheme =
> MatrixFree<dim, double>::AdditionalData::none;
> additional_data.mapping_update_flags =
> (update_values | update_gradients | update_JxW_values |
> update_quadrature_points);
> additional_data.mapping_update_flags_boundary_faces =
> (update_values | 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_matrix.initialize(system_mf_storage);
> system_mf_storage->initialize_dof_vector(alpha_solution,1);
> system_mf_storage->initialize_dof_vector(system_rhs,1);
> }
>
> On Thursday, June 22, 2023 at 11:59:31 AM UTC-6 Sean Johnson wrote:
>
>> I can further confirm the problem is with what is called in step-37
>> mg_matrices. Since in setup the mg_matrices is aware of both the CG
>> dof_handler and DG dof_handler it leads to it making a larger size.
>> However, the diagonal_inverse vector is setup with the CG dof_handler which
>> in turn does not match the mg_matrices size when initializing in the
>> smoother when making the preconditioner.
>>
>> Currently, for system_matrix and mg_matrices I have simply only changed
>> the reinit() functions to account for a std::vector of DoFHandlers and a
>> std::vector of AffineConstraints. In addition, this leads to needing to
>> form a std::vector of MGConstrainedDoFs for the mg_matrices. I thought this
>> was needed so that I could use the values of the DG variables in the
>> computation of the system_matrix and mg_matrices. Is there a different way
>> I should be doing this?
>>
>> On Wednesday, June 21, 2023 at 12:29:11 PM UTC-6 Sean Johnson wrote:
>>
>>> I am trying to make a matrix free solver like in step-37; however, my
>>> equations involve variables that are solved with a DG DoFHandler. So for
>>> solving I changed within LaplaceProblem<dim>::setup_system() the
>>> following functions to account for multiple DoFHandler's:
>>> system_mf_storage->reinit()
>>> mg_mf_storage_level->reinit()
>>> which led to me having to edit mg_matrices[level
>>> <https://www.dealii.org/current/doxygen/deal.II/grid__out_8cc.html#a9082f945c1d289684d0bcd51ee08e11e>].initialize()
>>>
>>> to account for multiple mg_constrained_dofs.
>>>
>>> All of that was fine until I got to LaplaceProblem<dim>::solve(). I cant
>>> edit either of the following functions because they don't allow for
>>> multiple DoFHandler's:
>>> MGTransferMatrixFree<dim, float>
>>> <https://www.dealii.org/current/doxygen/deal.II/classMGTransferMatrixFree.html>
>>>
>>> mg_transfer(mg_constrained_dofs);
>>> mg_transfer.build(dof_handler);
>>>
>>> I am unsure if that will cause problems but the line that for sure does
>>> cause an error and for the code not to finish is initializing the
>>> mg_smoother:
>>> mg_smoother.initialize(mg_matrices, smoother_data);
>>>
>>> This always gives me the error:
>>> An error occurred in line <2684> of file
>>> </home/sjohnson/AE_beginnings/dealii/include/deal.II/lac/precondition.h> in
>>> function
>>> void
>>> dealii::internal::PreconditionChebyshevImplementation::initialize_preconditioner(const
>>>
>>> MatrixType&, std::shared_ptr<dealii::DiagonalMatrix<VectorType> >&) [with
>>> MatrixType = Brill_Evolution::alpha_Operator<2, 4, float>; VectorType =
>>> dealii::LinearAlgebra::distributed::Vector<float>]
>>> The violated condition was:
>>> preconditioner->m() == 0
>>> Additional information:
>>> Preconditioner appears to be initialized but not sized correctly
>>>
>>> I believe this is because the mg_matrices has been initialized with
>>> multiple dof_handlers but I am unsure how to address this problem and have
>>> not been able to find a matrixfree example with multiple dof_handlers.
>>>
>>> Thanks,
>>> Sean Johnson
>>>
>>
--
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/45123ba0-0f26-4810-8659-1329eb9c7d75n%40googlegroups.com.
