Matrix-free DG with h-adaptive mesh refinement

[M. Bänsch]

Dear community,

I am currently trying to implement h-adaptive mesh refinement into step-67 (by adding a refine_grid() functionality). As references I looked into step-31, step-33, the previous question for step-59 and the distributed::SolutionTransfer documentation. 
I tried to adapt the approach in the SolutionTransfer documentation so that it matches what I see in the references. 


Unfortunately, I get the following error:

-----------------
The violated condition was:
    (data_range.size() % bytes_per_entry == 0)
Additional information:
    This exception -- which is used in many places in the library -- usually indicates that some condition which the author of the code thought must be satisfied at a certain point in an algorithm, is not fulfilled. An example would be that the first part of an algorithm sorts elements of an array in ascending order, and a second part of the algorithm later encounters an element that is not larger than the previous one.

There is usually not very much you can do if you encounter such an exception since it indicates an error in deal.II, not in your own program. Try to come up with the smallest possible program that still demonstrates the error and contact the deal.II mailing lists with it to obtain help.
-----------------

The previous question basically achieved what I am trying to implement but the additional functionality that comes with step-59's setup_system() call takes care of all the interpolation etc.

I have the feeling that I misunderstood something in the documentation (for example replacing the interpolated_solution part. Would someone be able to tell me if I made a mistake with adapting the documentation to what I have. 
Is there also some mistake on my part since I need another structure for LinearAlgebra::distributed::Vector (instead of the TrilinosWrappers vectors from the documentation)?

Any help would be much appreciated.

Best regards,
M. Bänsch

This is the code snippet:

template <int dim>
  void EulerProblem<dim>::refine_grid()
  {
...
some refinement criteria
...

// see distributed::SolutionTransfer documentation for details
  parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<Number>>     soltrans(dof_handler);
  triangulation.prepare_coarsening_and_refinement();
  soltrans.prepare_for_coarsening_and_refinement(solution);
  triangulation.execute_coarsening_and_refinement();

  dof_handler.distribute_dofs(fe);
  IndexSet locally_owned_dofs, locally_relevant_dofs;
  locally_owned_dofs = dof_handler.locally_owned_dofs();
  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);

  LinearAlgebra::distributed::Vector<Number> solution;
  solution.reinit(locally_owned_dofs, MPI_COMM_WORLD);

  LinearAlgebra::distributed::Vector<Number> old_solution;
  old_solution.reinit(locally_owned_dofs,
                    locally_relevant_dofs,
                    MPI_COMM_WORLD);
  old_solution = solution;
  soltrans.prepare_for_coarsening_and_refinement(old_solution);

  solution.reinit(locally_owned_dofs, MPI_COMM_WORLD);
soltrans.interpolate(solution);

}

[Wolfgang Bangerth]

On 8/8/22 02:44, M. Bänsch wrote:
> *** Caution: EXTERNAL Sender ***

>
> Dear community,
>
> I am currently trying to implement h-adaptive mesh refinement into step-67 (by
> adding a refine_grid() functionality). As references I looked into step-31

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fstep_31.html&data=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7C748797de48544ab12a7f08da791a3eb8%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637955450882849163%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=6qPfXFxiq0ozBmOtX8MYP5BrnmuvzYY7iBJLeM%2Bv%2BYg%3D&reserved=0>,
> step-33
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fstep_33.html&data=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7C748797de48544ab12a7f08da791a3eb8%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637955450882849163%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=XCKIzscMt9PFmPe5aUYV3%2FNSuOBsPCCYk5QgcoR1TIw%3D&reserved=0>,

> the previous question for step-59

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fg%2Fdealii%2Fc%2FokNWVIHwcTA%2Fm%2FqhfflrXkBgAJ&data=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7C748797de48544ab12a7f08da791a3eb8%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637955450882849163%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=7m%2BsPB5lHeunKIo%2BmzlIdYkXBRpfVehRIge3MijHYEU%3D&reserved=0>

> and the distributed::SolutionTransfer documentation.
> I tried to adapt the approach in the SolutionTransfer documentation

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fclassparallel_1_1distributed_1_1SolutionTransfer.html&data=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7C748797de48544ab12a7f08da791a3eb8%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637955450882849163%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=SpPh8sR2M%2Fd6YD1gngnxbZFCHzcRKbj5spcwGm1RUn4%3D&reserved=0>

It's hard to say what the reason is without having a small test case that
shows the issue. Since the issue happens during solution transfer, it is
almost certainly not something related to whether or not you build solutions.

As is generically the case, it is easiest to find the problem if you reduce
the program that produces the error to something as minimal as possible --
just throw out everything that can be thrown out. Do you think you can come up
with such a program?

Separately (and almost certainly unrelated): Did you mean to reinit() the
'solution' variable twice here?
```

LinearAlgebra::distributed::Vector<Number> old_solution;
old_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
MPI_COMM_WORLD);
old_solution = solution;
soltrans.prepare_for_coarsening_and_refinement(old_solution);

solution.reinit(locally_owned_dofs, MPI_COMM_WORLD);
soltrans.interpolate(solution);

```

Best
W.


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

[M. Bänsch]

On Wednesday, August 10, 2022 at 1:06:58 AM UTC+2 Wolfgang Bangerth wrote:

It's hard to say what the reason is without having a small test case that
shows the issue. Since the issue happens during solution transfer, it is
almost certainly not something related to whether or not you build solutions.

As is generically the case, it is easiest to find the problem if you reduce
the program that produces the error to something as minimal as possible --
just throw out everything that can be thrown out. Do you think you can come up
with such a program?

 

Thank you very much for your input! Your suggestion regarding the minimal working example is much appreciated. I trimmed down a code snippet even more and the current version is attached. This example is roughly 340 lines long and the relevant part with refine_grid() start at line 204. I included a dummy refinement criterion. 

 

 

Separately (and almost certainly unrelated): Did you mean to reinit() the
'solution' variable twice here?
```
LinearAlgebra::distributed::Vector<Number> old_solution;
old_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
MPI_COMM_WORLD);
old_solution = solution;
soltrans.prepare_for_coarsening_and_refinement(old_solution);

solution.reinit(locally_owned_dofs, MPI_COMM_WORLD);
soltrans.interpolate(solution);
```

For this section I just followed the documentation for the distributed SolutionTransfer but I was also a little bit confused about this part, to be honest. I am currently playing around with this part to see what the different outcomes can be.

Best regards,
M. Bänsch 

[Wolfgang Bangerth]

On 8/13/22 11:20, M. Bänsch wrote:
> *** Caution: EXTERNAL Sender ***
>
>

> On Wednesday, August 10, 2022 at 1:06:58 AM UTC+2 Wolfgang Bangerth wrote:
>
>
> It's hard to say what the reason is without having a small test case that
> shows the issue. Since the issue happens during solution transfer, it is
> almost certainly not something related to whether or not you build solutions.
>
> As is generically the case, it is easiest to find the problem if you reduce
> the program that produces the error to something as minimal as possible --
> just throw out everything that can be thrown out. Do you think you can
> come up
> with such a program?
>
> Thank you very much for your input! Your suggestion regarding the minimal
> working example is much appreciated. I trimmed down a code snippet even more
> and the current version is attached. This example is roughly 340 lines long
> and the relevant part with refine_grid() start at line 204. I included a dummy
> refinement criterion.

I got it down to just over 100 lines and put the result into this bug report:
https://github.com/dealii/dealii/issues/14189

Unfortunately, I don't have the time right now to debug this, but at least
it's recorded in our bug database. I also tagged it to make sure we fix this
before the release. I recognize that that isn't helpful to you right now, so
if you want to debug this yourself: you might want to read through the file
where the error happens and try to understand what the exception here is
trying to assert:

(data_range.size() % bytes_per_entry == 0)

This is in a function called unpack_dof_values(), and I imagine it is checking
whether some data it gets has been packed in the way that this function
expects, but that that check fails for some reason. Once you understand *what*
the assertion checks, you can start to ask *why* the assertion fails in your case.

[Marc Fehling]

Hello!

You are calling prepare_for_coarsening_and_refinement() twice which is not allowed. Thanks for reporting this on the mailing list -- we will improve the error message that you encountered.

I spotted a few issues in the order of events of your program. Tutorials step-32, step-42, step-48, step-70 demonstrate the use of parallel::distributed::SolutionTransfer class. With the help of these, I am sure you will get your program to run!

Marc