underworldcode / underworld3

https://underworldcode.github.io/underworld3/
18 stars 10 forks source link

Unexpected kernel crash with essential boundary conditions #221

Open gthyagi opened 1 month ago

gthyagi commented 1 month ago

Hi @lmoresi @julesghub @knepley ,

Issue: It is currently not possible to implement a essential boundary condition using a mesh variable that is not available in the Stokes object. In the example below, v_uw is the primary mesh variable available in the Stokes object, while v_ana is the mesh variable used just to store the top boundary velocities.

Issue lies here: https://github.com/underworldcode/underworld3/blob/d0b6f3b71ca9dea282571e7ad470c7de6b90925f/src/underworld3/cython/petsc_generic_snes_solvers.pyx#L2720C9-L2750C1

Here the block of code to replicate the essential boundary issue.

import underworld3 as uw
from underworld3.systems import Stokes
import sympy

mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1)

# mesh variables
v_uw = uw.discretisation.MeshVariable('V_u', mesh, mesh.data.shape[1], degree=2)
p_uw = uw.discretisation.MeshVariable('P_u', mesh, 1, degree=1, continuous=True)
v_ana = uw.discretisation.MeshVariable('V_ana', mesh, mesh.data.shape[1], degree=2)

with mesh.access(v_ana):
    v_ana.data[:,0] = 1.0

# Create Stokes object
stokes = Stokes(mesh, velocityField=v_uw, pressureField=p_uw, solver_name="stokes_sph")
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.viscosity = 1
stokes.saddle_preconditioner = 1.0/1.0

# +
# gravity
gravity_fn = -1.0 * mesh.CoordinateSystem.unit_e_1

# bodyforce term
stokes.bodyforce = 1.0*gravity_fn # 0.0 * unit_rvec
# -

# bc's
stokes.add_essential_bc(v_ana.sym, "Top") # v_ana.sym
stokes.add_essential_bc((sympy.oo,0.0), "Bottom")
stokes.add_essential_bc((0.0,sympy.oo), "Left")
stokes.add_essential_bc((0.0,sympy.oo), "Right")

stokes.solve(verbose=True, debug=False)
knepley commented 1 month ago
  1. The only fields available in kernels are either solution fields (U, P) or auxiliary fields, but I do not see V_ana marked as auxiliary.
  2. I do not understand the BC on V_ana. Enforcing a boundary condition on a non-solution field does not make sense to me.
julesghub commented 1 month ago

Thanks for the upload @gthyagi. For now if you replace stokes.add_essential_bc(v_ana.sym, "Top") with stokes.add_essential_bc((sympy.oo, 1), "Top") does that work?

gthyagi commented 1 month ago

yes, stokes.add_essential_bc((sympy.oo, 1), "Top") works.

gthyagi commented 1 month ago

@knepley, V_ana is used to store the velocity analytical solution for the benchmark problems and to impose the tangential component of boundary velocities on the inner/outer surface. Currently, essential boundary conditions can accept either constant numerical values or sympy expressions for each component of velocity. However, if I provide a mesh variable like V_ana.sym or its components (i.e., [V_ana.sym[0], V_ana.sym[1]]) to set the boundary values, an error occurs.

gthyagi commented 1 week ago

@lmoresi @knepley I'm trying to debug this issue and have included the details of the error below. Could you suggest any possible fixes based on the provided information?

Here are the UW3 lines: https://github.com/gthyagi/underworld3/blob/2f985eaf20a1f1e4bd91d98f5ecb5e9ad8541ca6/src/underworld3/cython/petsc_generic_snes_solvers.pyx#L2721-L2750

Error happens in the line where DMPlexSNESComputeBoundaryFEM() called.

I am using petsc 3.21.1.

SNES solve - picard = 0
SNES post-solve - bcs
Processing field pressure
Processing field pressure
Obtained subvector for pressure
Obtained subvector for pressure
Local and sub-dm vectors obtained for pressure
Local and sub-dm vectors obtained for pressure
Field - sgvec size: 142, lvec size: 79
Field - sgvec size: 142, lvec size: 76
Field - sgvec size: 142, lvec size: 79
Field - sgvec size: 142, lvec size: 76
Field - sgvec size: 142, lvec size: 76
Field - sgvec size: 142, lvec size: 79
Inserting boundary conditions for pressure
Inserting boundary conditions for pressure
Field pressure processing completed.
Processing field velocity
Field pressure processing completed.
Processing field velocity
Obtained subvector for velocity
Obtained subvector for velocity
Local and sub-dm vectors obtained for velocity
Local and sub-dm vectors obtained for velocity
Field - sgvec size: 947, lvec size: 538
Field - sgvec size: 947, lvec size: 562
Field - sgvec size: 947, lvec size: 562
Field - sgvec size: 947, lvec size: 538
Field - sgvec size: 947, lvec size: 538
Field - sgvec size: 947, lvec size: 562
Inserting boundary conditions for velocity
Inserting boundary conditions for velocity
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
[1]PETSC ERROR: ------------------------------------------------------------------------
[1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[1]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[1]PETSC ERROR: ---------------------  Stack Frames ------------------------------------
[1]PETSC ERROR: The line numbers in the error traceback are not always exact.
[1]PETSC ERROR: #1 DMProjectPoint_Field_Private() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/impls/plex/plexproject.c:298
[1]PETSC ERROR: #2 DMProjectPoint_Private() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/impls/plex/plexproject.c:485
[1]PETSC ERROR: #3 DMProjectLocal_Generic_Plex() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/impls/plex/plexproject.c:970
[1]PETSC ERROR: #4 DMProjectFieldLabelLocal_Plex() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/impls/plex/plexproject.c:1053
[1]PETSC ERROR: #5 DMProjectFieldLabelLocal() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/interface/dm.c:8368
[1]PETSC ERROR: #6 DMPlexInsertBoundaryValuesEssentialField() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/impls/plex/plexfem.c:953
[1]PETSC ERROR: #7 DMPlexInsertBoundaryValues_Plex() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/impls/plex/plexfem.c:1173
[1]PETSC ERROR: #8 DMPlexInsertBoundaryValues() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/dm/impls/plex/plexfem.c:1274
[1]PETSC ERROR: #9 DMPlexSNESComputeBoundaryFEM() at /Users/tgol0006/manual_install_pkg/venv_uw3_git_gthyagi_5924/petsc/src/snes/utils/dmplexsnes.c:469
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 0 on node MU00155712X exited on signal 11 (Segmentation fault: 11).
--------------------------------------------------------------------------