GalerkinToolkit.jl

GalerkinToolkit.jl copied to clipboard

As can be seen below, and based on our discussion on 2024-04-02, the sequential two level mesh for the periodic puzzle piece appears to have issues at the interface. The visualize_mesh helper function, and # TODO: comments in example.jl commit 1fa5b349c814dbcb02d4d6feccac3600cde9d17e indicate areas of the code that produce the preliminary diagnostic information we discussed today. The top of example.jl describes the naming conventions, but for reading the below files it's worth noting final_mesh vtk files are named as final_mesh_unit_cell_mesh_<info about unit cell mesh>_coarse_cell_mesh_<info about coarse cell mesh>.

The image below is from the output/final_mesh_unit_cell_mesh_3D_periodic_gmsh_puzzlepiece_geometry_triangular_refcell_coarse_cell_mesh_3D_periodic_glk_box_geometry_quad_2x2x2_refcell.vtu.

It is also worth inspecting output/coarse_cell_id_final_mesh_unit_cell_mesh_3D_periodic_gmsh_box_geometry_triangular_refcell_coarse_cell_mesh_3D_periodic_glk_box_geometry_quad_4x4x4_refcell.vtu since the coarse cell ids appear to be labeled correctly.

Perhaps I found the issue. See this comment in the code: https://github.com/fverdugo/GalerkinToolkit.jl/blob/eed905c038c9356653f8c51db40696340d14e500/src/geometry.jl#L3644

We cannot assume a single level of indirection in 3D. We need 2 levels of indirection. In general, we need D-1 levels of indirection.

Let's try this and check if it fixes the issue:

for _ in 1:(D-1)
    for fine_node in 1:n_fine_nodes
        master_node = fine_node_to_master_node[fine_node]
        if fine_node == master_node
            continue  
        end
        master_master_node = fine_node_to_master_node[master_node]
        fine_node_to_master_node[fine_node] = master_master_node 
    end
end

Commit a73aaa9cce2eb9f618a729e761a763c4304635a5, particularly src/tmp.jl is an attempt to handle the master-slave problem introduced by 3D unit cell geometries. The image below is a screenshot of outputs/debug_3D_periodic_puzzlepiece.vtu, which can be reproduced by running include("example.jl") from the project main directory. The bracketed edges represent problematic edges that, when inspecting mesh nodes, reveal that both edges are treated as masters, though it is expected that at least one edge should be a slave. Note also the axial orientation in the bottom left corner.

Good news!

@jfdev001 the correct implementation is in a function I created for debugging here:

https://github.com/fverdugo/GalerkinToolkit.jl/blob/d6f91b3374421a56fc4981fca7963e3cfa9b2ced/src/geometry.jl#L3999

can you please move the fix to the two_level_mesh function and delete the new function?

After this, also check the function in parallel. If this work, we are ready for clean up, a pull request, and start running experiments!

I can explain you which was the problem in the next meeting.

@fverdugo I have updated the src/geometry.jl function in commit 63d431814ab66d801514b4b71a1157f5c468a6d7 to use the corrected version you provided. Corresponding tests in example.jl all pass; however, I think we need to still have tests to verify node and cell ownership, see some logic considerations in

https://github.com/fverdugo/GalerkinToolkit.jl/blob/63d431814ab66d801514b4b71a1157f5c468a6d7/example.jl#L1286

and

https://github.com/fverdugo/GalerkinToolkit.jl/blob/63d431814ab66d801514b4b71a1157f5c468a6d7/example.jl#L1294

I think these tests are still necessary because using paraview does not necessarily make this clear in the 3D case, and hovering over the points with the inspection tool is not tractable for these complicated meshes.

Lastly, I have added a small solver test in

https://github.com/fverdugo/GalerkinToolkit.jl/blob/63d431814ab66d801514b4b71a1157f5c468a6d7/extensions/GalerkinToolkitExamples/test/example002_defs.jl#L75

which is called accordingly in

https://github.com/fverdugo/GalerkinToolkit.jl/blob/63d431814ab66d801514b4b71a1157f5c468a6d7/extensions/GalerkinToolkitExamples/test/mpi_array/example002_driver_np_4.jl#L7

but the solver fails to converge for a 2D periodic puzzle piece here

https://github.com/fverdugo/GalerkinToolkit.jl/blob/63d431814ab66d801514b4b71a1157f5c468a6d7/extensions/GalerkinToolkitExamples/test/example002_defs.jl#L137

function label_boundary_faces!(mesh::PMesh;physical_name="boundary")
    D = num_dims(mesh)
    face_parts = face_partition(mesh,D-1)
    v = pfill(1,face_parts)
    assemble!(v) |> wait
    map(partition(mesh),partition(v),face_parts) do mymesh, myv, myfaces
        topo = topology(mymesh)
        face_to_cells = face_incidence(topo,d,D)
        local_to_owner_face = local_to_owner(myfaces)
        part = part_id(myfaces)
        ids = findall(1:length(myfaces)) do face
            myv[face] == 1 && local_to_owner_face[face] == part && length(face_to_cells) == 1
        end
        physical_faces(mymesh,D-1)[physical_name] = ids
    end
    mesh
end

@fverdugo 5bdd15afc12aa98c7c90adf2141f377d643fdd6e incorporates the new function

https://github.com/fverdugo/GalerkinToolkit.jl/blob/5bdd15afc12aa98c7c90adf2141f377d643fdd6e/src/geometry.jl#L3450

and makes some modifications to correspond with the logic in

https://github.com/fverdugo/GalerkinToolkit.jl/blob/5bdd15afc12aa98c7c90adf2141f377d643fdd6e/src/geometry.jl#L2246

To debug, I wrote

https://github.com/fverdugo/GalerkinToolkit.jl/blob/5bdd15afc12aa98c7c90adf2141f377d643fdd6e/extensions/GalerkinToolkitExamples/example.jl#L42

and while that tolerance test now passes, the paraview visualization seems to disagree with the logic we discussed yesterday (i.e., the "boundary" tag should correspond to parts of the mesh that are not on the interface between partitions). To reproduce the below visualization, do the following:

cd extensions/GalerkinToolkitExamples
julia
julia> ]
(@<julia version>) pkg> activate .
(GalerkinToolkitExamples) pkg> <backspace>
julia> include("example.jl")
julia> TMP.test_solver_periodic_2D_puzzle_piece_pmesh()

and see GalerkinToolkitExamples/output/puzzle_piece_2D_np_4_test.pvtu.

Per a0ad02d692d28393b22cd6ad32bddd4bec8d5deb, all solver tests for 2D square unit cells, 3D box unit cells, 2D puzzle piece unit cell, and 3D puzzle piece unit cell in both sequential and parallel pass. Upon visualizing the boundary in GalerkinToolkitExamples/output/puzzle_piece_3D_np_4_test.pvtu, reproducible following

cd extensions/GalerkinToolkitExamples
julia
julia> ]
(@<julia version>) pkg> activate .
(GalerkinToolkitExamples) pkg> <backspace>
julia> include("example.jl")

the boundary visualization looks off, see below:

Per a previous comment,

I think we need to still have tests to verify node and cell ownership

because this is critical to parallel GalerkinToolkit functionality, though I have not looked too far into this yet.

Hi @jfdev001 Try to use the implementation in my last comment. (I introduced some changes just after our meeting).

https://github.com/fverdugo/GalerkinToolkit.jl/issues/65#issuecomment-2058598819

Commit 9dc9880858c269d3eb36ff5d34273318f81c2847 was the first time I incorporated your function from https://github.com/fverdugo/GalerkinToolkit.jl/issues/65#issuecomment-2058598819; however, the tests for the 2D pmesh puzzle piece fail for GalerkinToolkitExamples/example.jl and TMP.test_solver_periodic_2D_puzzle_piece_pmesh(). I saw that the boundaries were not actually be annotated with your original function implementation, as shown below in GalerkinToolkitExamples/output/puzzle_piece_2D_np_4_test.pvtu

and so I updated your function as in https://github.com/fverdugo/GalerkinToolkit.jl/issues/65#issuecomment-2060923981

Now I remember that we just defined a "dummy" face partition in the parallel two_level_mesh, right?

I think that there is nothing wrong with my label_boundary_faces! funciton. The problem is that the face partition in the mesh is wrong.

@fverdugo As of our discussion on 2024-04-23, the use of dummy face partitions for faces in 0:D-1 leads to a failure of the label boundary faces algorithm. Attempts to fix this problem are on a new branch called broken-label-boundary-faces, the latest commit for which is fb45ed1a43a461d708763cef84c446da9d473b3f. The branch parallel_two_level_mesh uses a modified, but incorrect, label_boundary_faces! function which allows all solver tests to pass in GalerkinToolkitExamples/example.jl. I will clean up the parallel_two_level_mesh branch, move its tests into the appropriate files, and then submit a pull request noting that label_boundary_faces! in the parallel case is incorrectly implemented to the main branch so that I can move on to setting up experiments.