Help setting periodic boundary in box

[bknight1]

Trying to get a periodic boundary using the following:

_mesh = uw3.meshing.StructuredQuadBox(minCoords=(0, 0), maxCoords=(1, 1), elementRes=(36, 36) )

plex = uw3.discretisation._from_gmsh(".meshes/uw_structuredQuadBox_minC(0, 0)_maxC(1, 1).msh",
                useMultipleTags=True,
                useRegions=True,)

uw3.cython.petsc_discretisation.petsc_dm_set_periodicity(
     plex[1], [1.0, 0.0], [0.0, 0.0], [1.0, 0.0])

### match boundary labels from plex to the boundaries class
from enum import Enum

class boundaries_2D(Enum):
    Top = 12
    Right = 13
    Bottom = 11
    Left = 14

mesh = uw3.discretisation.Mesh(plex[1], useMultipleTags=True, useRegions=True, 
                              boundaries=boundaries_2D, coordinate_system_type=uw3.coordinates.CoordinateSystemType.CARTESIAN)

However, the solution does not appear to have a periodic boundary condition. Am I setting it up correctly?

[lmoresi]

Have a look at the example: "Examples-StokesFlow/Ex_WIP_Stokes_Periodic.py"

It has the gmsh stuff and the dmplex setting-up stuff that you need. I'm not sure how we would handle particle advection or particle-based semi-Lagrange advection across that wraparound with the current setup but that may not matter for your case.

[bknight1]

Thanks Louis. I resorted to using PETSc for now, which seems to be working when no velocity is imposed (Ice benchmark).

### Use PETSc to create periodic mesh

from petsc4py import PETSc
plex = PETSc.DMPlex().createBoxMesh(
    lower=[minX, minY],  # Lower corner of the mesh (x, y)
    upper=[maxX, maxY],  # Upper corner of the mesh (x, y)
    faces=[60, 24],  # Number of faces along each axis
    simplex=False,  # Use quadrilateral (False) or triangular (True) elements
    periodic=[True, False]  # Periodicity in each dimension (x, y)
)

# Define markers for each boundary
markers_dict = {"Bottom": 1,
                "Right": 2,
                "Top": 3,
                "Left": 4}

# Define tolerance for boundary detection
tol = 1e-20

# Get the mesh coordinates and reshape them for easier access
arr = plex.getCoordinates().array.reshape(-1, 2)

# Loop through the markers and assign labels
for key, value in markers_dict.items():
    plex.createLabel(key)
    label = plex.getLabel(key)

    if key == "Bottom":
        # Find points where y is close to the minimum y value (bottom boundary)
        indices = [i for i, coord in enumerate(arr) if abs(coord[1] - arr[:, 1].min()) < tol]
    elif key == "Top":
        # Find points where y is close to the maximum y value (top boundary)
        indices = [i for i, coord in enumerate(arr) if abs(coord[1] - arr[:, 1].max()) < tol]
    elif key == "Left":
        # Find points where x is close to the minimum x value (left boundary)
        indices = [i for i, coord in enumerate(arr) if abs(coord[0] - arr[:, 0].min()) < tol]
    elif key == "Right":
        # Find points where x is close to the maximum x value (right boundary)
        indices = [i for i, coord in enumerate(arr) if abs(coord[0] - arr[:, 0].max()) < tol]

    # print(len(indices))

    # Convert indices to PETSc IS
    indexSet = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF)

    # Insert indices into the label
    label.insertIS(indexSet, value)

    # Cleanup
    indexSet.destroy()

# # Print labels for verification
# for key in markers_dict:
print(f"Label '{key}' has points:", len(plex.getStratumIS(key, markers_dict[key]).getIndices()) if plex.getStratumIS(key, markers_dict[key]) else "None")

# ### match boundary labels from plex to the boundaries class
from enum import Enum

class boundaries_2D(Enum):
    Top = 3
    Right = 2
    Bottom = 1
    Left = 4

mesh = uw.discretisation.Mesh(plex, useMultipleTags=True, useRegions=True, 
                              boundaries=boundaries_2D, coordinate_system_type=uw.coordinates.CoordinateSystemType.CARTESIAN)

However, seems to be an issue when using velocity BC with period conditions, will be investigating further

[lmoresi]

Is it still the case that PETSc cannot export to VTK for our visualisation module when periodic ?

[knepley]

If Plex is making the mesh, -dm_plex_periodic_cut is supposed to mark a line in the mesh that we can cut to unroll it, using a DMLabel. Then the viz part cuts it open when generating output. We can make this work for an arbitrary mesh if we mark the cut line I think.

On Thu, Aug 29, 2024 at 8:05 AM Louis Moresi @.***> wrote:

Is it still the case that PETSc cannot export to VTK for our visualisation module when periodic ?

— Reply to this email directly, view it on GitHub https://github.com/underworldcode/underworld3/issues/236#issuecomment-2317457377, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAEORCPVRFFSSGL3QMVVBCTZT4E7LAVCNFSM6AAAAABNHMW3CGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJXGQ2TOMZXG4 . You are receiving this because you are subscribed to this thread.Message ID: @.***>

-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ http://www.cse.buffalo.edu/~knepley/

[bknight1]

I had an issue with higher order velocity fields when using the integration points with matplotlib, so no saving as a vtk. The velocity vectors from the periodic boundary vectors were close to the centre of the plot, plotting over another set of coordinates.

when v_deg = 3 (continuous) and p_deg = 2 (discontinuous):

when v_deg = 2 (continuous) and p_deg = 1 (discontinuous):

when v_deg = 1 (continuous) and p_deg = 0 (discontinuous):

When trying to save the mesh created by DMPlex as a vtk I get the following error:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: VTK output with localized coordinates not yet supported
[0]PETSC ERROR: WARNING! There are unused option(s) set! Could be the program crashed before usage or a spelling mistake, etc!

[lmoresi]

Yes, that PETSC ERROR was also my experience. However, for the purposes of visualisation, the vtk mesh can be constructed on the fly (see what I am doing for higher order fields). We dump the fields with swarm-like information for exactly this kind of reason. I don't really care about the detailed connectivity in most cases, just the values, locations.

The same issue with the erroneous location of some values - it's clearly just a storage / post-processing issue so we can work around it [For now ...]

[bknight1]

Closing this as the initial issue was solved