step-37 multiple DoFHandler's smoother problems

[Sean Johnson]

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].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> 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

[Sean Johnson]

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?

[Sean Johnson]

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

[Sean Johnson]

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

[Sean Johnson]

I keep seeing other questions getting at least one response. Is there something else that I should provide in addition to the question or that is missing or is there somewhere this question has already been addressed that I cannot find? I am really all ears

[Wolfgang Bangerth]

On 6/27/23 09:53, Sean Johnson wrote:
> I keep seeing other questions getting at least one response. Is there
> something else that I should provide in addition to the question or that is
> missing or is there somewhere this question has already been addressed that I
> cannot find? I am really all ears

Sean, it's all about who has the expertise to answer questions and who has the
time. I don't know anything about matrix-free stuff, so your questions require
someone else to have the time and expertise.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Sean Johnson]

Okay. Thank you.

[Peter Munch]

I'll try to take a look at the issue tomorrow. Could you attach the full code. I have a feeling what might be an issue. 

Peter

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/dab3cc77-bd64-445a-82f5-da55b00a812en%40googlegroups.com.

[Sean Johnson]

Yes of course. I can get one with minimal changes made to step-37. Full code do you want the entirety of the changed step-37 code or just the LaplaceProblem class or just the setup_system() function?

[Sean Johnson]

Here is step-37 with minimal changes. All that was changed was in the creation of a DG DoFHandler and the fe_system and constrains and what not to go along with it. Then of the functions only the LapalceProblem setup_system() function was altered and nothing else. It was just altered to also be reinit with the dof_handler_DG.

Obviously the LaplaceOperator in my actual code actually uses variables that depend on the DG DoFHandler but that was left out of this example.

Thank you again so much for your time and for all you guys do for this library.

Thanks,

Sean Johnson

[Peter Munch]

OK. The problem are these lines: https://github.com/dealii/dealii/blob/8d68b3192ed1c8b928dfce03ca1d615201859904/include/deal.II/matrix_free/operators.h#L1225-L1235. The number of rows are the sum of the number of DoFs of both attached DoFHandlers...

Personally, I would not use these matrix-free operators but create the operator myself around a MatrixFree object. Just like in step-75. Would that be possible?

Out of curiosity, what is your goal?

Best,
PM

[Sean Johnson]

I haven't done step-75 but I will look into it today and take a shot at it. There is nothing special about the way I did it besides that step-37 came earlier so I had already done it and was familiar with that method. Thank you for your time and suggestion. I'll give an update probably in a day or so after learning and trying to do step-75.

Big picture goal is we are doing a gravity waves simulation where they implode and form a black hole. Numerically, it has a few variables that have equations for evolving forward in time and then like three that don't but do have elliptic equations that if we know all the variables that have time equations we can then solve them on each time step. Thus the need for a DG DoFHandler for the time equations and a normal CG DoFHandler for the elliptic equation variables.

Thanks again so much for your time,

Sean Johnson

[Sean Johnson]

Okay I think I have implemented what I could of step-75. I currently have deal.II compiled without Trilinos and am trying to keep it that way because it is compiled without Trilinos on the high powered computer we will do runs on later. That can change but would cause a headache for someone else so I am just trying to avoid it.

So I created the operator around a MatrixFree object. It allowed me the ability to control m() for the matrix. I still have a problem and I went back to my old version and this was a problem in both of them. Since in the reinit function of the MatrixFree object is given two DoFHandlers, I have a system matrix that expects a two component solution vector and rhs. The solution vector and rhs are scalars with just one component. I am handing both DoFHandlers over just because the computation of the matrix elements depend on some variables I have calculated elsewhere with another DoFHandler.

Thanks for your time so far. I am still actively working on how to address this. If you have any suggestions I am all ears.

I will also work on a minimal example that causes the error.

Thanks for your time again,

Sean Johnson

[Sean Johnson]

I believe all I had to do was upon initializing the MatrixFree operator I needed to include a vector that specified the "selected_row_blocks" which allows you to choose which component to use in the underlying matrix free object by including a vector with the index of the dof_handlers and constriants you choose to use. The n_components function I used still returned that the object had 2 components but the solving went through as if there were only one.

Thanks for your patience,

Sean Johnson