Matrix Free 1D Adaptive Mesh Problem

[Sean Johnson]

Hey,

It might be my fault but I made a simplest case pair of files attached.

In 1D for adaptive meshes, the matrix free operator is not iterating over the face/point that is between refinement levels.

I attached a code that basically just makes a mesh and refines the last cell. I use the dof_handler to print out all the coordinates of the faces. It does iterate over all the faces. However, the matrix free operator in its face iterator doesn't recognize the point between refinement levels as a face so I have it print out all the coordinates of faces and its obviously missing one.

The "mesh" is just a line from [0,16]. First it is divided into 4 cells. Then the final cell is further refined so there should be faces at: 0,4,8,12,14,16.

Again, feel free to correct me if I made a dumb mistake somewhere, but I did try the code in 2D and there were no problems there. The code is super simplistic and will not give you enough information about coordinates if you step up to 3D.

Thanks for any help you can give,

Sean Johnson

[Wolfgang Bangerth]

On 11/5/24 13:53, Sean Johnson wrote:
>
> In 1D for adaptive meshes, the matrix free operator is not iterating
> over the face/point that is between refinement levels.
>
> I attached a code that basically just makes a mesh and refines the last
> cell. I use the dof_handler to print out all the coordinates of the
> faces. It does iterate over all the faces. However, the matrix free
> operator in its face iterator doesn't recognize the point between
> refinement levels as a face so I have it print out all the coordinates
> of faces and its obviously missing one.
>
> The "mesh" is just a line from [0,16]. First it is divided into 4 cells.
> Then the final cell is further refined so there should be faces at:
> 0,4,8,12,14,16.
>
> Again, feel free to correct me if I made a dumb mistake somewhere, but I
> did try the code in 2D and there were no problems there. The code is
> super simplistic and will not give you enough information about
> coordinates if you step up to 3D.

Sean:
I don't know much about the matrix-fee framework, and so can't say much
about whether it does what it's supposed to -- someone else will have to
chime in here.

But can you describe what the code you attach outputs, and how that
differs from what you *expect* it to output? It's often difficult for
those who didn't write the code to say what is expected, and how what
you get differs from the expectation.

Best
WB

[Sean Johnson]

Wolfgang,

Thanks for responding. Of course I can explain. So as stated my code creates a mesh from 0 to 16 with 5 cells. First from 0 to 4, then 4 to 8, then 8 to 12, then 12 to 14, and finally 14 to 16.

The code first outputs all the coordinate of all the face values just by iterating over them using the dof_handler cell/face iterator. So it prints:

############# Printing X values of center of faces ################
Face at x_val: 0
Face at x_val: 4
Face at x_val: 4
Face at x_val: 8
Face at x_val: 8
Face at x_val: 12
Face at x_val: 12
Face at x_val: 14
Face at x_val: 14
Face at x_val: 16

Then, it creates a matrix free operator and goes through applying the matrix free operator to a vector of zeros. The matrix free operator is built to do nothing for the cell and boundary face loops. However, for the interior face loops it is just supposed to print the coordinates of the quadrature points (since the face is just a point in 1D it should print the coordinate of the face). What it prints is:

#######################################################
Matrix Free Face Iteration
#######################################################
x_val: 4
x_val: 8
x_val: 14
x_val: 14

Then just to verify that the matrix free operator didn't change anything and I didn't have stray code changing anything I repeat the first part of the code that uses the dof_handler and it prints the exact same thing.

So what I expect to be different is this second part with the matrix free operator. I would expect it to print:

#######################################################
Matrix Free Face Iteration
#######################################################
x_val: 4
x_val: 8
x_val: 12
x_val: 14

Or basically to have "x_val: 12" anywhere in this section because 12 is the coordinate of an interior face of the mesh and is not currently included in the output of the code.. This only happens at the face between two cells of different refinement levels and only in 1D from my testing.

Thanks for your time,

Sean Johnson

[Sean Johnson]

Hey,

Sorry I am still digging through trying to find out why this happens.

I have confirmed that the matrix free operator thinks there are 2 batches of interior faces with the first one containing a pair of interior faces and the second just containing one face. There are 4 pairs of interior faces so this still implies that deal.ii in matrix free 1D with an adaptive mesh misses the face between refinement levels.

I was trying to read the code of how the matrix free operator is initialized but it uses a function called internal_reinit() and I can't find the actual code defining this function anywhere in matrix_free.h

Any help would be greatly appreciated.

Thanks,

Sean Johnson

[Sean Johnson]

I have found what I believe to be the cause of the problem. Its in the dealii file face_setup_internal.h which is inside the matrix_free folder.

To account for adaptive meshes, when creating faces the files checks if the neighbor has any children to know if it might be further refined. If so, it then passes through the children of the face and flips a flag on them. The problem is in 1D the cell does have children; however, the faces return 0 for the number of children so no flags get flipped. This is all around line 220 of the file. In addition, the create_face function in the same file then also needs a slight tweak in the if statement around line 994 that has the condition "cell->level() > neighbor->level()". The else statement (which will be where a normal 1D adaptive mesh face will go to) uses a function that is not intended for 1D meshes. I just changed it to an else if statement that required the dim to be larger than 1.

I am in the process of testing if this works in an actual problem being solved. Previously, matrix free adaptive meshes just would never evaluate the integral on faces between different refinement levels which obviously causes large errors. The solution mentioned above fixed the issue in my simplest version of the problem script that I provided in the first message. The script just printed which faces were visited by the the matrix_free operator and which faces did exist in the dof_handler.

I'll update if it doesn't fully resolve my original issue.

Sean Johnson

[Sean Johnson]

I misspoke. The first problem isn't line 220 it is around 816.

It didn't fully fix the problem as one of the values was populated with the wrong value but I would bet this has to do with my renumbering of the Degrees of Freedom.

[Wolfgang Bangerth]

Sean,
I am not familiar with the matrix-free framework, but I wanted to say that if
you figure out how to make this work, we would certainly very much like it if
you could submit your fixes as a patch to deal.II! It sounds like you're doing
a good job working your way through the code and understanding what needs to
be fixed -- that is work that would be useful for others as well!

Best
Wolfgang


On 12/9/24 15:45, Sean Johnson wrote:
> *** Caution: EXTERNAL Sender ***

> --
> The deal.II project is located at http://www.dealii.org/ <https://
> nam10.safelinks.protection.outlook.com/?
> url=http%3A%2F%2Fwww.dealii.org%2F&data=05%7C02%7CWolfgang.Bangerth%40colostate.edu%7C98c10b8299c94e8832dc08dd18a341f1%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C638693811665413564%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=oAjy4MZnaqM%2FCKOxcpuTnaJzwWDiSaICKcldkqNcicI%3D&reserved=0>
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?
> hl=en <https://nam10.safelinks.protection.outlook.com/?
> url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=05%7C02%7CWolfgang.Bangerth%40colostate.edu%7C98c10b8299c94e8832dc08dd18a341f1%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C638693811665439568%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=hgFwYYkb54lIVcVqvK98Zl%2BlTQnfSopvDU01J6f2n60%3D&reserved=0>
> ---
> 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
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion visit https://groups.google.com/d/msgid/dealii/
> e28280cc-0893-481d-bafa-f5e9c7fac310n%40googlegroups.com <https://
> nam10.safelinks.protection.outlook.com/?
> url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2Fe28280cc-0893-481d-
> bafa-
> f5e9c7fac310n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=05%7C02%7CWolfgang.Bangerth%40colostate.edu%7C98c10b8299c94e8832dc08dd18a341f1%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C638693811665457012%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=txgCCeeoByPKMB6JcnE3PzH2fgibyycWVWb4SUXmvDU%3D&reserved=0>.

[Sean Johnson]

Wolfgang,

Thanks! I did end up getting it all working and I have now submitted a pull request to dealii.

Thanks again for all the hard work you all have put into this!

Best,

Sean Johnson