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

Hybrid Jacobian integration gives incorrect values for cell matrix #103

Closed baagaard-usgs closed 5 years ago

baagaard-usgs commented 5 years ago

The fault kernels appear to be called correctly and the values are added into the cell Jacobian matrix, but the entries are canceling each other out, resulting in a matrix full of zeros.

These are constraint equations, so the entries should not cancel each other out.

baagaard-usgs commented 5 years ago

The indexing issue that cause the values to cancel is gone.

The hybrid residual integration seems to be working. When I set the initial conditions to the solution (displacement and Lagrange multiplier values), the residual is zero and SNES correctly determines that the solution has converged.

In the case where I do not give an initial condition (start with zeros as the initial guess for the solution), SNES fails to converge.

Test case (no initial conditions)

PyLith branch: baagaard/feature-integrator-interface

cd tests_auto/linearelasticity/faults-2d
make export-data
pylith twoblocks.cfg twoblocks_quad.cfg

Solver output

0 TS dt 0.01 time 0.
    0 SNES Function norm 4.377975178855e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 5
      0 KSP Residual norm 3.626411640193e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 5
      1 KSP Residual norm 2.465396306957e-15
    Linear solve converged due to CONVERGED_ATOL iterations 1
        Line search: Using full step: fnorm 4.377975178855e-03 gnorm 2.213901656402e-03
    1 SNES Function norm 2.213901656402e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 5
      0 KSP Residual norm 1.992377892726e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 6
      1 KSP Residual norm 2.172185169766e-16
    Linear solve converged due to CONVERGED_ATOL iterations 1
        Line search: Using full step: fnorm 2.213901656402e-03 gnorm 1.417708027154e-03
    2 SNES Function norm 1.417708027154e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 6
      0 KSP Residual norm 2.342192632652e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 6
      1 KSP Residual norm 6.958095615048e-15
    Linear solve converged due to CONVERGED_ATOL iterations 1
        Line search: Using full step: fnorm 1.417708027154e-03 gnorm 1.210521518432e-03
    3 SNES Function norm 1.210521518432e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 6
      0 KSP Residual norm 3.468741067459e-03
      Linear fieldsplit_lagrange_multiplier_fault_ solve converged due to CONVERGED_RTOL iterations 6
      1 KSP Residual norm 5.943303775305e-15
    Linear solve converged due to CONVERGED_ATOL iterations 1
        Line search: gnorm after quadratic fit 1.264150339899e-03
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.224969223260e-03 lambda=1.6347578580665778e-01
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.215399071397e-03 lambda=6.6312608731712658e-02
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.212345954162e-03 lambda=2.7027050071017997e-02
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.211241142447e-03 lambda=1.1064686786300080e-02
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210811954619e-03 lambda=4.5360813382502672e-03
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210639891941e-03 lambda=1.8608104543732505e-03
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210569960599e-03 lambda=7.6354093712843339e-04
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210541375818e-03 lambda=3.1333448801297283e-04
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210529663983e-03 lambda=1.2858865172849930e-04
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210524860706e-03 lambda=5.2772143874668305e-05
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210522889990e-03 lambda=2.1657582018543169e-05
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210522081301e-03 lambda=8.8882545783450739e-06
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521749430e-03 lambda=3.6477371050846466e-06
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521613233e-03 lambda=1.4970309990018221e-06
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521557338e-03 lambda=6.1438152102725782e-07
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521534399e-03 lambda=2.5214219699556780e-07
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521524985e-03 lambda=1.0347916978532104e-07
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521521121e-03 lambda=4.2467856898716172e-08
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521519535e-03 lambda=1.7428810853962870e-08
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518885e-03 lambda=7.1527850037541680e-09
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518618e-03 lambda=2.9355033984872975e-09
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518508e-03 lambda=1.2047308710621793e-09
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518463e-03 lambda=4.9442169136336799e-10
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518445e-03 lambda=2.0291073769094289e-10
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518437e-03 lambda=8.3274693660356366e-11
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518434e-03 lambda=3.4175988425609224e-11
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518433e-03 lambda=1.4025809998483626e-11
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518432e-03 lambda=5.7562265886543245e-12
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518432e-03 lambda=2.3623741425849645e-12
        Line search: Cubic step no good, shrinking lambda, current gnorm 1.210521518432e-03 lambda=9.6953132199175330e-13
        Line search: unable to find good step length! After 30 tries
        Line search: fnorm=1.2105215184317814e-03, gnorm=1.2105215184318428e-03, ynorm=3.4687410674585670e-03, minlambda=9.9999999999999998e-13, lambda=9.6953132199175330e-13, initial slope=-1.4653623465863740e-06
  Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH iterations 3

Note: I have snes_max_it = 4, but this is a linear problem so it should converge in 1 iteration.

baagaard-usgs commented 5 years ago

If I let the solution be written by setting

snes_error_if_not_converged = false
ts_error_if_step_fails = false

I see that the displacements are zero except on the -x and +x boundaries where I have Dirichlet BC. I am not sure if this is problem is due to incorrect solver settings or an incorrect Jacobian or other.

baagaard-usgs commented 5 years ago

FIxed. More testing needed.

baagaard-usgs commented 5 years ago

Matt's fix for buried fault edges messed up faults without buried edges.

I just pushed a fix to baagaard/reorganize-examples so that the constraint on fault edges is not created when they don't exist.

However, the Jacobian is messed up and I get "KSPSolve has not converged due to Nan or Inf norm" when running any problem with a through-going fault (i.e., no buried edges). Setting dm_plex_print_fem = 2 shows garbage values in the fault Jacobian. The residual appears to be ok.

Test case (pull/build baagaard/reorganize-examples)

cd $BUILD_DIR/tests_auto/linearelasticity/faults-2d
make export-data
pylith twoblocks.cfg twoblocks_quad.cfg --petsc.dm_plex_print_fem=2
knepley commented 5 years ago

@baagaard I just ran this. At line 218 of FaultCohesiveKin.cc, which is

err = DMGetLabel(dm, getBuriedEdgesMarkerLabel(), &buriedLabel);PYLITH_CHECK_ERROR(err);

I get that buriedLabel is missing. Ah, which should be true if there is no buried edge. Fixing now.

knepley commented 5 years ago

@baagaard I made a branch knepley/fix-fault-noburied for this from baagaard/regorganized-examples. What is the workflow now?

baagaard-usgs commented 5 years ago

So you think you have a fix, but we need to test it, right? If you push your branch, Charles and I can run tests and if the buried edges / non-buried edges are fixed, we merge to baagaard/reorganize-examples.

knepley commented 5 years ago

The branch is pushed.

baagaard-usgs commented 5 years ago

@knepley Your branch is from June 13 and is out-of-date with the current reorganize-examples. I believe I already have a better fix for what you found in 2ada3a9875d75e9e71978aaf4f2492d450e236d6.

I am running more tests as faults with and without buried edges may be working again.

baagaard-usgs commented 5 years ago

PyLith branch: baagaard/reorganize-examples

I get the same behavior as I did on Jun 28 (See https://github.com/geodynamics/pylith/issues/103#issuecomment-506845887).

PyLith branch: knepley/fix-fault-noburied

I get the same behavior as before.

On my Linux machine the garbage values include very large numbers so the solve fails, but on my Mac the garbage values are very small, so the solve proceeds.

cd $BUILD_DIR/tests_auto/linearelasticity/faults-2d
make export-data
pylith twoblocks.cfg twoblocks_quad.cfg --petsc.dm_plex_print_fem=2
--snip--
Cell 16 Element Hybrid Jacobian
  |         0         0         0         0         0         0         0         0   0.66667  2.7401e-233   0.33333  4.9631e-234 |
  |         0         0         0         0         0         0         0         0  3.0479e+287   0.66667  3.0479e+287   0.33333 |
  |         0         0         0         0         0         0         0         0   0.33333  2.7401e-233   0.66667  4.9631e-234 |
  |         0         0         0         0         0         0         0         0 -3.0479e+287   0.33333 -3.0479e+287   0.66667 |
  |         0         0         0         0         0         0         0         0  -0.66667  8.9403e-25  -0.33333 -8.9403e-25 |
  |         0         0         0         0         0         0         0         0  2.4756e+287  2.8965e+226  2.4756e+287 -2.8965e+226 |
  |         0         0         0         0         0         0         0         0  -0.33333  8.9403e-25  -0.66667 -8.9403e-25 |
  |         0         0         0         0         0         0         0         0 -2.4756e+287 -2.8965e+226 -2.4756e+287  2.8965e+226 |
  |   0.66667  2.7401e-233   0.33333  4.9631e-234  -0.66667 -8.0381e-102  -0.33333 -1.1314e-266         0         0         0         0 |
  |  5.7237e+286 -2.8965e+226  5.7237e+286  2.8965e+226  4.7362e-218 -1.0564e+106  4.745e-218  1.0564e+106         0         0         0         0 |
  |   0.33333  2.7401e-233   0.66667  4.9631e-234  -0.33333 -8.0381e-102  -0.66667 -1.1314e-266         0         0         0         0 |
  | -3.0479e+287   0.33333 -3.0479e+287   0.66667 -1.596e-204  -0.33333 -1.596e-204  -0.66667         0         0         0         0 |
  Indices: 6 7 0 1 32 33 30 31 42 43 40 41

Mac output

Cell 16 Element Hybrid Jacobian
  |         0         0         0         0         0         0         0         0   0.66667  4.2384e-319   0.33333 -1.1348e-314 |
  |         0         0         0         0         0         0         0         0  4.1905e-319   0.66667  4.1905e-319   0.33333 |
  |         0         0         0         0         0         0         0         0   0.33333  4.2384e-319   0.66667 -1.1348e-314 |
  |         0         0         0         0         0         0         0         0 -1.1348e-314   0.33333 -1.1348e-314   0.66667 |
  |         0         0         0         0         0         0         0         0  -0.66667  3.3102e-322  -0.33333 -1.1347e-314 |
  |         0         0         0         0         0         0         0         0 -1.5934e-320  -0.66667 -1.5934e-320  -0.33333 |
  |         0         0         0         0         0         0         0         0  -0.33333  3.3102e-322  -0.66667 -1.1347e-314 |
  |         0         0         0         0         0         0         0         0 -1.1347e-314  -0.33333 -1.1347e-314  -0.66667 |
  |   0.66667  4.2384e-319   0.33333 -1.1348e-314  -0.66667  4.5503e-321  -0.33333 -1.1347e-314         0         0         0         0 |
  |  6.1306e-319   0.66667  6.1306e-319   0.33333  1.7808e-319  -0.66667  1.7808e-319  -0.33333         0         0         0         0 |
  |   0.33333  4.2384e-319   0.66667 -1.1348e-314  -0.33333  4.5503e-321  -0.66667 -1.1347e-314         0         0         0         0 |
  | -1.1348e-314   0.33333 -1.1348e-314   0.66667 -1.1347e-314  -0.33333 -1.1347e-314  -0.66667         0         0         0         0 |
  Indices: 6 7 0 1 32 33 30 31 42 43 40 41
knepley commented 5 years ago

Okay, I think I have a fix. Pushed to knepley/pylith

baagaard-usgs commented 5 years ago

I must have been running something old yesterday, because now I get the same random values. Valgrind shows a number of errors. See attached logs (I have include cases both with and without the dm_plex_pring_fem=2).

cd tests_auto/linearelasticity/faults-2d
make export-data
valgrind --trace-children=yes --suppressions=$HOME/src/cig/pylith/share/valgrind-python.supp pylith twoblocks.cfg twoblocks_quad.cfg --petsc.dm_plex_print_fem=2 >& valgrind.log

valgrind.log valgrind_noprint.log

knepley commented 5 years ago

Okay, I found one thing but it is not the problem. There is one Constraint that is being applied over a boundary mesh, but the type that gets passed in is ESSENTIAL_FIELD, which causes things to mess up (I think). I believe it should have type ESSENTIAL_BD_FIELD.

baagaard-usgs commented 5 years ago

All of our Dirichlet BC (feassemble/ConstraintSimple.cc, feassemble/ConstraintSpatialDB.cc, and feassemble/ConstraintUserFn.cc) use DM_BC_ESSENTIAL or DM_BC_ESSENTIAL_FIELD. Nowhere in PyLith do we use just ESSENTIAL_FIELD.

knepley commented 5 years ago

Is it possible we have used DM_NATURAL_FIELD? I have this stack trace:

==43305== Invalid read of size 8 ==43305== at 0x8BF5372: PetscFEEvaluateFieldJets_Internal (fe.c:1737) ==43305== by 0x8AEA1B7: DMProjectPoint_Field_Private (plexproject.c:221) ==43305== by 0x8AE760E: DMProjectPoint_Private (plexproject.c:368) ==43305== by 0x8AE35D0: DMProjectLocal_Generic_Plex (plexproject.c:701) ==43305== by 0x8AE5BF7: DMProjectFieldLabelLocal_Plex (plexproject.c:828) ==43305== by 0x8C5E185: DMProjectFieldLabelLocal (dm.c:7807) ==43305== by 0x8A56DC1: DMPlexInsertBoundaryValuesEssentialField (plexfem.c:799) ==43305== by 0x76695BD: pylith::feassemble::ConstraintSpatialDB::setSolution(pylith::topology::Field, double) (in /PETSc3/cig/lib/libpylith.0.dylib) ==43305== by 0x7739F89: pylith::problems::Problem::setSolutionLocal(double, _p_Vec, _p_Vec) (in /PETSc3/cig/lib/libpylith.0.dylib) ==43305== by 0x773BB36: pylith::problems::Problem::computeLHSResidual(_p_Vec, double, double, _p_Vec, _p_Vec) (in /PETSc3/cig/lib/libpylith.0.dylib) ==43305== by 0x774E332: pylith::problems::TimeDependent::computeLHSResidual(_p_TS, double, _p_Vec, _p_Vec, _p_Vec, void*) (in /PETSc3/cig/lib/libpylith.0.dylib) ==43305== by 0x93B3D0F: TSComputeIFunction (ts.c:891) ==43305== Address 0x1050dcc48 is 4 bytes after a block of size 1,892 alloc'd ==43305== at 0x47E1: malloc (vg_replace_malloc.c:300) ==43305== by 0x7A8EAD0: PetscMallocAlign (mal.c:62) ==43305== by 0x7A91FD4: PetscTrMallocDefault (mtr.c:190) ==43305== by 0x8C18A00: DMGetWorkArray (dm.c:1383) ==43305== by 0x8BEA16F: PetscFEGetTabulation (fe.c:893) ==43305== by 0x8AE1EE9: DMProjectLocal_Generic_Plex (plexproject.c:607) ==43305== by 0x8AE5BF7: DMProjectFieldLabelLocal_Plex (plexproject.c:828) ==43305== by 0x8C5E185: DMProjectFieldLabelLocal (dm.c:7807)

baagaard-usgs commented 5 years ago

We never use DM_NATURAL_FIELD. We call the boundary residual functions explicitly.

I don't understand what you are asking about the stack trace?

knepley commented 5 years ago

This point in the stack trace:

DMProjectPoint_Private (plexproject.c:368)

is pointing here in the code

362 switch (type) { 363 case DM_BC_ESSENTIAL: 364 case DM_BC_NATURAL: 365 ierr = DMProjectPoint_Func_Private(dm, ds, dmIn, dsIn, time, fegeom, &fvgeom, isFE, sp, (PetscErrorCode ()(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar , void )) funcs, ctxs, values);CHKERRQ(ierr);break; 366 case DM_BC_ESSENTIAL_FIELD: 367 case DM_BC_NATURAL_FIELD: 368 ierr = DMProjectPoint_Field_Private(dm, ds, dmIn, encIn, dsIn, dmAux, encAux, dsAux, time, localU, localA, fegeom, sp, p, basisTab, basisDerTab, basisTabAux, basisDerTabAux, (void ()(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) funcs, ctxs, values);CHKERRQ(ierr);break;

So it certainly looks like one of those values is being used.

knepley commented 5 years ago

@baagaard Okay, in setSolution() we call DMPlexInsertBoundaryValuesEssentialField() which uses DMProjectFieldLabelLocal() which uses DM_BC_ESSENTIAL_FIELD by default, even though in this case we really want DM_BC_ESSENTIAL_BD_FIELD. I guess its time to rethink this part of the interface since now lots of different integrations are possible.

baagaard-usgs commented 5 years ago

@knepley Are you saying we need to change PetscDSAddBoundary(prob, DM_BC_ESSENTIAL_FIELD, ...) to PetscDSAddBoundary(prob, DM_BC_ESSENTIAL_BD_FIELD, ...) and then make sure DMProjectFieldLabelLocal() uses this flag rather than assuming DM_BC_ESSENTIAL_FIELD?

When I use PetscDSAddBoundary(prob, DM_BC_ESSENTIAL_BD_FIELD, ...) I see DMProjecPoint_Private() is still using type=DM_BC_ESSENTIAL_FIELD.

knepley commented 5 years ago

Yes, I think that sounds like the right solution.

knepley commented 5 years ago

Finally did this and pushed to knepley/pylith

knepley commented 5 years ago

@baagaard I just pushed knepley/pylith and made a branch knepley/fix-bc-kernels from baagaard/reorganize-examples (I could not branch from master). I am now valgrind clean on faults-2d/twoblocks_quad. Let me know if it works for you

baagaard-usgs commented 5 years ago

I am getting errors when I try to run other problems.

PyLith branch: knepley/fix-bc-kernels

tests_auto/linearelasticity/nofaults-2d
pylith axialdisp.cfg axialdisp_tri.cfg
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: No label named depth was found
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.12-124-g0eccbdd1fc  GIT Date: 2019-10-04 12:31:07 -0400
[0]PETSC ERROR: /tools/common/cig/gcc-7.3.0/bin/mpinemesis on a arch-gcc-7.3.0_debug named arling by brad Fri Oct  4 12:28:25 2019
[0]PETSC ERROR: Configure options --PETSC_ARCH=arch-gcc-7.3.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 --CC=mpicc --CXX=mpicxx --FC= --FC= LDFLAGS= --with-blaslapack-lib="-L/usr/lib/atlas-base -llapack_atlas -llapack -latlas -lblas" --with-lgrind=0 --with-c2html=0 --with-chaco=1 --download-chaco=1 --download-parmetis=1 --download-metis=1 --with-ml=1 --download-ml=1 --download-superlu --with-hdf5=1 --with-zlib=1 --with-hdf5-dir=/tools/common/hdf5-1.10.2/gcc-7.3.0 --download-triangle --with-fc=0
[0]PETSC ERROR: #1 DMPlexGetHeightStratum() line 3661 in /tools/common/petsc-dev/src/dm/impls/plex/plex.c
[0]PETSC ERROR: #2 DMFieldCreateDS() line 1081 in /tools/common/petsc-dev/src/dm/field/impls/ds/dmfieldds.c
[0]PETSC ERROR: #3 DMCreateCoordinateField_Plex() line 3703 in /tools/common/petsc-dev/src/dm/impls/plex/plex.c
[0]PETSC ERROR: #4 DMGetCoordinateField() line 5743 in /tools/common/petsc-dev/src/dm/interface/dm.c
[0]PETSC ERROR: #5 DMProjectLocal_Generic_Plex() line 587 in /tools/common/petsc-dev/src/dm/impls/plex/plexproject.c
[0]PETSC ERROR: #6 DMProjectFunctionLocal_Plex() line 793 in /tools/common/petsc-dev/src/dm/impls/plex/plexproject.c
[0]PETSC ERROR: #7 DMProjectFunctionLocal() line 7729 in /tools/common/petsc-dev/src/dm/interface/dm.c
[0]PETSC ERROR: #8 void pylith::topology::FieldQuery::queryDB()() line 214 in /home/brad/src/cig/pylith/libsrc/pylith/topology/FieldQuery.cc
Fatal error. Calling MPI_Abort() to abort PyLith application.
Traceback (most recent call last):
  File "/tools/common/pylith-develop/gcc-7.3.0/lib/python2.7/site-packages/pylith/apps/PetscApplication.py", line 74, in onComputeNodes
    self.main(*args, **kwds)
  File "/tools/common/pylith-develop/gcc-7.3.0/lib/python2.7/site-packages/pylith/apps/PyLithApp.py", line 119, in main
    self.problem.initialize()
  File "/tools/common/pylith-develop/gcc-7.3.0/lib/python2.7/site-packages/pylith/problems/Problem.py", line 201, in initialize
    ModuleProblem.initialize(self)
  File "/tools/common/pylith-develop/gcc-7.3.0/lib/python2.7/site-packages/pylith/problems/problems.py", line 145, in initialize
    def initialize(self): return _problems.Problem_initialize(self)
RuntimeError: Error detected while in PETSc function.
knepley commented 5 years ago

Well, I think I have most of the bugs fixed from this. The newest knepley/pylith should be enough. However, right now I am getting some leaks. I cannot understand why I do not get a stack. I will keep looking...

baagaard-usgs commented 5 years ago

The unit tests in unittests/libtests/topology show a leak and you can run the valgrind leak check using make leakcheck.

==15357== Conditional jump or move depends on uninitialised value(s)
==15357==    at 0x6912EB3: CompressPoints_Private (plex.c:4228)
==15357==    by 0x6913E2F: DMPlexGetCompressedClosure (plex.c:4311)
==15357==    by 0x6915BA1: DMPlexVecGetClosure (plex.c:4498)
==15357==    by 0x1A0637: getClosure (CoordsVisitor.icc:120)
==15357==    by 0x1A0637: pylith::topology::TestReverseCuthillMcKee::testReorder() (TestReverseCuthillMcKee.cc:119)
==15357==    by 0x1A52EE: __invoke_impl<void, void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:73)
==15357==    by 0x1A52EE: __invoke<void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:95)
==15357==    by 0x1A52EE: __call<void, 0> (functional:467)
==15357==    by 0x1A52EE: operator()<> (functional:551)
==15357==    by 0x1A52EE: std::_Function_handler<void (), std::_Bind<void (pylith::topology::TestReverseCuthillMcKee::*(pylith::topology::TestReverseCuthillMcKee*))()> >::_M_invoke(std::_Any_data const&) (std_function.h:316)
==15357==    by 0x135A1E: std::function<void ()>::operator()() const (std_function.h:706)
==15357==    by 0x1AB971: CppUnit::TestCaller<pylith::topology::TestReverseCuthillMcKee>::runTest() (TestCaller.h:175)
==15357==    by 0x4E5E881: CppUnit::TestCaseMethodFunctor::operator()() const (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E54E62: CppUnit::DefaultProtector::protect(CppUnit::Functor const&, CppUnit::ProtectorContext const&) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E5B9E4: CppUnit::ProtectorChain::protect(CppUnit::Functor const&, CppUnit::ProtectorContext const&) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E6441B: CppUnit::TestResult::protect(CppUnit::Functor const&, CppUnit::Test*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E5E6BF: CppUnit::TestCase::run(CppUnit::TestResult*) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==  Uninitialised value was created by a heap allocation
==15357==    at 0x4C320A6: memalign (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==15357==    by 0x5A30109: PetscMallocAlign (mal.c:49)
==15357==    by 0x6C813C0: DMGetWorkArray (dm.c:1380)
==15357==    by 0x68FF47A: DMPlexGetTransitiveClosure (plex.c:2344)
==15357==    by 0x52BFBF6: pylith::topology::MeshOps::isSimplexMesh(pylith::topology::Mesh const&) (MeshOps.cc:225)
==15357==    by 0x52BD3C3: pylith::topology::Mesh::dmMesh(_p_DM*, char const*) (Mesh.cc:98)
==15357==    by 0x5158F61: pylith::faults::TopologyOps::create(pylith::topology::Mesh*, pylith::topology::Mesh const&, _p_DMLabel*, int) (TopologyOps.cc:205)
==15357==    by 0x51397D8: pylith::faults::FaultCohesive::adjustTopology(pylith::topology::Mesh*) (FaultCohesive.cc:203)
==15357==    by 0x19E059: pylith::topology::TestReverseCuthillMcKee::_initialize() (TestReverseCuthillMcKee.cc:218)
==15357==    by 0x19E4FB: pylith::topology::TestReverseCuthillMcKee::testReorder() (TestReverseCuthillMcKee.cc:65)
==15357==    by 0x1A52EE: __invoke_impl<void, void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:73)
==15357==    by 0x1A52EE: __invoke<void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:95)
==15357==    by 0x1A52EE: __call<void, 0> (functional:467)
==15357==    by 0x1A52EE: operator()<> (functional:551)
==15357==    by 0x1A52EE: std::_Function_handler<void (), std::_Bind<void (pylith::topology::TestReverseCuthillMcKee::*(pylith::topology::TestReverseCuthillMcKee*))()> >::_M_invoke(std::_Any_data const&) (std_function.h:316)
==15357==    by 0x135A1E: std::function<void ()>::operator()() const (std_function.h:706)
==15357== 
==15357== Conditional jump or move depends on uninitialised value(s)
==15357==    at 0x6912EBB: CompressPoints_Private (plex.c:4228)
==15357==    by 0x6913E2F: DMPlexGetCompressedClosure (plex.c:4311)
==15357==    by 0x6915BA1: DMPlexVecGetClosure (plex.c:4498)
==15357==    by 0x1A0637: getClosure (CoordsVisitor.icc:120)
==15357==    by 0x1A0637: pylith::topology::TestReverseCuthillMcKee::testReorder() (TestReverseCuthillMcKee.cc:119)
==15357==    by 0x1A52EE: __invoke_impl<void, void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:73)
==15357==    by 0x1A52EE: __invoke<void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:95)
==15357==    by 0x1A52EE: __call<void, 0> (functional:467)
==15357==    by 0x1A52EE: operator()<> (functional:551)
==15357==    by 0x1A52EE: std::_Function_handler<void (), std::_Bind<void (pylith::topology::TestReverseCuthillMcKee::*(pylith::topology::TestReverseCuthillMcKee*))()> >::_M_invoke(std::_Any_data const&) (std_function.h:316)
==15357==    by 0x135A1E: std::function<void ()>::operator()() const (std_function.h:706)
==15357==    by 0x1AB971: CppUnit::TestCaller<pylith::topology::TestReverseCuthillMcKee>::runTest() (TestCaller.h:175)
==15357==    by 0x4E5E881: CppUnit::TestCaseMethodFunctor::operator()() const (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E54E62: CppUnit::DefaultProtector::protect(CppUnit::Functor const&, CppUnit::ProtectorContext const&) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E5B9E4: CppUnit::ProtectorChain::protect(CppUnit::Functor const&, CppUnit::ProtectorContext const&) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E6441B: CppUnit::TestResult::protect(CppUnit::Functor const&, CppUnit::Test*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==    by 0x4E5E6BF: CppUnit::TestCase::run(CppUnit::TestResult*) (in /usr/lib/x86_64-linux-gnu/libcppunit-1.14.so.0.0.0)
==15357==  Uninitialised value was created by a heap allocation
==15357==    at 0x4C320A6: memalign (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==15357==    by 0x5A30109: PetscMallocAlign (mal.c:49)
==15357==    by 0x6C813C0: DMGetWorkArray (dm.c:1380)
==15357==    by 0x68FF47A: DMPlexGetTransitiveClosure (plex.c:2344)
==15357==    by 0x52BFBF6: pylith::topology::MeshOps::isSimplexMesh(pylith::topology::Mesh const&) (MeshOps.cc:225)
==15357==    by 0x52BD3C3: pylith::topology::Mesh::dmMesh(_p_DM*, char const*) (Mesh.cc:98)
==15357==    by 0x5158F61: pylith::faults::TopologyOps::create(pylith::topology::Mesh*, pylith::topology::Mesh const&, _p_DMLabel*, int) (TopologyOps.cc:205)
==15357==    by 0x51397D8: pylith::faults::FaultCohesive::adjustTopology(pylith::topology::Mesh*) (FaultCohesive.cc:203)
==15357==    by 0x19E059: pylith::topology::TestReverseCuthillMcKee::_initialize() (TestReverseCuthillMcKee.cc:218)
==15357==    by 0x19E4FB: pylith::topology::TestReverseCuthillMcKee::testReorder() (TestReverseCuthillMcKee.cc:65)
==15357==    by 0x1A52EE: __invoke_impl<void, void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:73)
==15357==    by 0x1A52EE: __invoke<void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:95)
==15357==    by 0x1A52EE: __call<void, 0> (functional:467)
==15357==    by 0x1A52EE: operator()<> (functional:551)
==15357==    by 0x1A52EE: std::_Function_handler<void (), std::_Bind<void (pylith::topology::TestReverseCuthillMcKee::*(pylith::topology::TestReverseCuthillMcKee*))()> >::_M_invoke(std::_Any_data const&) (std_function.h:316)
==15357==    by 0x135A1E: std::function<void ()>::operator()() const (std_function.h:706)

==15357== 32,196 (5,120 direct, 27,076 indirect) bytes in 8 blocks are definitely lost in loss record 541 of 545
==15357==    at 0x4C320A6: memalign (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==15357==    by 0x5A30109: PetscMallocAlign (mal.c:49)
==15357==    by 0x5A31AC4: PetscMallocA (mal.c:422)
==15357==    by 0x66B876A: DMLabelCreate (dmlabel.c:35)
==15357==    by 0x66BEEA1: DMLabelDuplicate (dmlabel.c:526)
==15357==    by 0x66C8B37: DMLabelPermute (dmlabel.c:1328)
==15357==    by 0x6A3EC13: DMPlexPermute (plexreorder.c:197)
==15357==    by 0x52F03A0: pylith::topology::ReverseCuthillMcKee::reorder(pylith::topology::Mesh*) (ReverseCuthillMcKee.cc:39)
==15357==    by 0x19E716: pylith::topology::TestReverseCuthillMcKee::testReorder() (TestReverseCuthillMcKee.cc:74)
==15357==    by 0x1A52EE: __invoke_impl<void, void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:73)
==15357==    by 0x1A52EE: __invoke<void (pylith::topology::TestReverseCuthillMcKee::*&)(), pylith::topology::TestReverseCuthillMcKee*&> (invoke.h:95)
==15357==    by 0x1A52EE: __call<void, 0> (functional:467)
==15357==    by 0x1A52EE: operator()<> (functional:551)
==15357==    by 0x1A52EE: std::_Function_handler<void (), std::_Bind<void (pylith::topology::TestReverseCuthillMcKee::*(pylith::topology::TestReverseCuthillMcKee*))()> >::_M_invoke(std::_Any_data const&) (std_function.h:316)
==15357==    by 0x135A1E: std::function<void ()>::operator()() const (std_function.h:706)
==15357==    by 0x1AB971: CppUnit::TestCaller<pylith::topology::TestReverseCuthillMcKee>::runTest() (TestCaller.h:175)