search

geodynamics / pylith

PyLith is a finite element code for the solution of dynamic and quasi-static tectonic deformation problems.
Other
154 stars 98 forks source link

PETSc error integrating boundary with hex cells #101

Closed baagaard-usgs closed 5 years ago

baagaard-usgs commented 5 years ago

PETSc error when integrating residual for a Neumann BC on a hex mesh.

PyLith branch: baagaard/feature-tests-auto-3d-nofault

Test case

Problem with tet mesh passes test. Hex mesh works fine with Dirichlet BC.

cd tests_auto/linearelasticity/nofaults-3d
make export-data
python sheartraction_gendb.py
pylith sheartraction.cfg sheartraction_hex.cfg

Output

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Cone size 6 not yet supported

[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.11.2-833-g688b753ea8  GIT Date: 2019-05-30 17:29:41 -0400
[0]PETSC ERROR: /Volumes/Tools/unix/cig/clang-3.6.0/bin/mpinemesis on a arch-clang-3.6.0_debug named IGSWMHWGLTE0904 by baagaard Tue Jun  4 21:01:12 2019
[0]PETSC ERROR: Configure options --PETSC_ARCH=arch-clang-3.6.0_debug --with-debugging=1 --with-clanguage=c --with-mpi-compilers=1 --with-shared-libraries=1 --with-64-bit-points=1 --with-large-file-io=1 --with-lgrind=0 --d
ownload-chaco=1 --download-parmetis=1 --download-metis=1 --download-triangle --with-blas-lib=/System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Versions/Current/libBLAS.dylib --with-lapack-lib=/Sys
tem/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Versions/Current/libLAPACK.dylib --download-ml=1 --download-superlu=1 --with-fc=0 --with-hdf5=1 --with-hdf5-include=/Volumes/Tools/unix/hdf5-1.10.5/cl
ang-3.6.0/include --with-hdf5-lib=/Volumes/Tools/unix/hdf5-1.10.5/clang-3.6.0/lib/libhdf5.dylib --with-zlib=1
[0]PETSC ERROR: #1 DMFieldComputeFaceData_DS() line 943 in /Volumes/Tools/unix/petsc-dev/src/dm/field/impls/ds/dmfieldds.c
[0]PETSC ERROR: #2 DMFieldCreateFEGeom() line 534 in /Volumes/Tools/unix/petsc-dev/src/dm/field/interface/dmfield.c
[0]PETSC ERROR: #3 DMSNESGetFEGeom() line 3506 in /Volumes/Tools/unix/petsc-dev/src/dm/impls/plex/plexfem.c
[0]PETSC ERROR: #4 DMPlexComputeBdResidual_Single_Internal() line 1118 in /Volumes/Tools/unix/petsc-dev/src/snes/utils/dmplexsnes.c
[0]PETSC ERROR: #5 DMPlexComputeBdResidualSingle() line 1206 in /Volumes/Tools/unix/petsc-dev/src/snes/utils/dmplexsnes.c
[0]PETSC ERROR: #6 static void pylith::feassemble::_IntegratorBoundary::computeResidual(pylith::topology::Field *, const pylith::feassemble::IntegratorBoundary *, const std::vector<pylith::feassemble::IntegratorBoundary::R
esidualKernels> &, const PylithReal, const PylithReal, const pylith::topology::Field &, const pylith::topology::Field &)() line 318 in /Volumes/Users/baagaard/src/cig/pylith/libsrc/pylith/feassemble/IntegratorBoundary.cc
baagaard-usgs commented 5 years ago

Update: Incorrect solution when using Neumann BC with hex cells.

PyLith branch: baagaard/reorganize-examples

Test case

Simulation with tet mesh passes test. Hex fails for both Dirichlet BC and Neumann BC.

cd tests_auto/linearelasticity/nofaults-3d
make export-data
nemesis TestShearTraction.py
knepley commented 5 years ago

@baagaard-usgs I just ran this and the Dirichlet tests for Hex also fail, which suggest it might be something other than the surface integration. Do you see that?

baagaard-usgs commented 5 years ago

In the shear traction example, we have both Dirichlet and Neumann BC. The Dirichlet BC do not constrain all three components, just the x and y DOF. This is why the error in the integration results in the z DOF being wrong. It still looks like the error is in the hex surface integration. The corresponding test in 2-D passes. That is, the 2-D test cases (nofaults-2d) all pass.

baagaard-usgs commented 5 years ago

@knepley says the face quadrature for hex cells is wrong. Neumann BCs work for tets and Matt says his face integration tests all use tets

baagaard-usgs commented 5 years ago

Need to remove all noncell info from material-id label after refining in RefineUniform.cc.

knepley commented 5 years ago

@baagaard I just tried to run this again to see the error, and I see no error if I run in the branch where I fixed refinement. Can you try running there to see if this indeed solves the issue?

baagaard-usgs commented 5 years ago

@knepley I still get the same error. Are you sure you ran the 3-D test in nofaults-3d? and not the 2-D test?

baagaard-usgs commented 5 years ago

Fixed in https://github.com/geodynamics/pylith/commit/ce362f517ac3adee7b673fdcdc2e34c00bb7bb31.