Refactor Method of Manufactured Solution tests

[baagaard-usgs]

PETSc branch: knepley/pylith PyLith branch: baagaard/feature-mms-testing

Use DMPlex routines for testing the residual and Jacobian via the Method of Manufactured Solutions.

@knepley needs to refactor his DMSNESCheckFromOptions_Internal() into pieces and return the norms to allow a clean implementation with CppUnit.

See DMTSCheckFromOptions(), DMSNESCheckFromOptions_Internal(), and nes/examples/tutorials/ex17.c.

Use functions for solution and boundary conditions via PetscDSSetExactSolution() and PetscDSAddBoundary().

[baagaard-usgs]

I have setup a sample MMS test case involving uniform strain for isotropic linear elasticity in 2D with plane strain. I currently call DMSNESCheckFromOptions_Internal() as a proxy for calling the refactored routines. See MMSTest.cc and TestIsotropicLinearElasticity2D_UniformStrain.cc.

I get the following error:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Argument aliasing not permitted
[0]PETSC ERROR: x, y, and z must be different vectors
-- snip --
[0]PETSC ERROR: #1 VecAXPBYPCZ() line 700 in /tools/common/petsc-dev/src/vec/vec/interface/rvector.c
[0]PETSC ERROR: #2 SNESTSFormFunction_Theta() line 644 in /tools/common/petsc-dev/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: #3 SNESTSFormFunction() line 4735 in /tools/common/petsc-dev/src/ts/interface/ts.c
[0]PETSC ERROR: #4 SNESComputeFunction() line 2406 in /tools/common/petsc-dev/src/snes/interface/snes.c
[0]PETSC ERROR: #5 DMSNESCheckFromOptions_Internal() line 2513 in /tools/common/petsc-dev/src/snes/utils/dmplexsnes.c

Reproduce

PETSc branch: knepley/pylith PyLith branch: baagaard/feature-mms-testing

cd mmstets/elasticity
make check

[baagaard-usgs]

Neither DMSNESCheckResidual() nor DMSNESCheckJacobian() set the solution at the beginning of the check. Refactoring DMSNESCheckFromOption() put setting the solution only in the DMSNESCheckDiscretization(). Because the residual and Jacobian checks require the solution to be set, I think it should be included there as well.

[knepley]

I am fixing this now. Does this sounds alright to you:

I will compute the exact solution if the solution function is in the DS. Otherwise, I will just take the values in the 'u' input

Thanks,

 Matt

On Mon, Apr 22, 2019 at 1:36 PM Brad Aagaard notifications@github.com wrote:

Neither DMSNESCheckResidual() nor DMSNESCheckJacobian() set the solution at the beginning of the check. Refactoring DMSNESCheckFromOption() put setting the solution only in the DMSNESCheckDiscretization(). Because the residual and Jacobian checks require the solution to be set, I think it should be included there as well.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/geodynamics/pylith/issues/72#issuecomment-485487708, or mute the thread https://github.com/notifications/unsubscribe-auth/AAEORCIYH6PNU6HNQLLJSG3PRXZRLANCNFSM4HDOTKGQ .

-- 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/

[baagaard-usgs]

Yes this seems reasonable.

I also want CppUnit to be managing the error reporting for these messages, so I want these functions to return some measure of the error (L2 diff, etc) that I can use to trigger a test pass/fail criterion.

[knepley]

On Mon, Apr 22, 2019 at 2:59 PM Brad Aagaard notifications@github.com wrote:

Yes this seems reasonable.

I also want CppUnit to be managing the error reporting for these messages, so I want these functions to return some measure of the error (L2 diff, etc) that I can use to trigger a test pass/fail criterion.

Okay, I have made both of these changes, and pushed to knepley/pylith

Thanks,

Matt

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/geodynamics/pylith/issues/72#issuecomment-485514260, or mute the thread https://github.com/notifications/unsubscribe-auth/AAEORCP5P6IBBDN4WNCNGFTPRYDKTANCNFSM4HDOTKGQ .

-- 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/

[baagaard-usgs]

MMS test fails for P3 and Q3.

PETSc branch: knepley/pylith PyLith branch: knepley/feature-petsc-fe

make install -C libsrc/pylith/
make check -C mmstests/elasticity

Test case: Uniform strain in two cell tri mesh with P3

pylith::mmstests::TestIsotropicLinearElasticity2D_UniformStrain_TriP3

Mesh

2 ---- 3
| \    |
|   \  |
|    \ |
0 ---- 1

cell 0: v0 v1 v2
cell 1: v1 v3 v2

Layout of solution

Subfields:
  Subfield displacement, index: 0, components: displacement_x displacement_y, scale: 1000, basisOrder: 3, quadOrder: 3, dimension: -1
  dimensionalize flag: 0
DM Object: 1 MPI processes
  type: plex
DM_0x55db714f1d20_83 in 2 dimensions:
  0-cells: 4
  1-cells: 5
  2-cells: 2
Labels:
  boundary: 1 strata with value/size (1 (9))
  material-id: 1 strata with value/size (24 (2))
  depth: 3 strata with value/size (0 (4), 1 (5), 2 (2))
Field displacement:
  adjacency FEM
PetscSection Object: 1 MPI processes
  type not yet set
1 fields
  field 0 with 2 components
Process 0:
  (   0) dim  2 offset   0
  (   1) dim  2 offset   2
  (   2) dim  2 offset   4 constrained 0 1
  (   3) dim  2 offset   6 constrained 0 1
  (   4) dim  2 offset   8 constrained 0 1
  (   5) dim  2 offset  10 constrained 0 1
  (   6) dim  4 offset  12 constrained 0 1
  (   7) dim  4 offset  16 constrained 0 1
  (   8) dim  4 offset  20 constrained 0 1
  (   9) dim  4 offset  24 constrained 0 1
  (  10) dim  4 offset  28 constrained 0 1
  PetscSectionSym Object: 1 MPI processes
    type: label
    Label 'depth'
    Symmetry for stratum value 0 (0 dofs per point): no symmetries
    Symmetry for stratum value 1 (4 dofs per point):
      Orientation range: [-2, 2)
    Symmetry for stratum value 2 (2 dofs per point):
      Orientation range: [-3, 3)
    Symmetry for stratum value -1 (0 dofs per point): no symmetries

Solution local vector has a size of 28. Solution global vector has a size of 14.

Shouldn't the edges on the boundary have 4 constraints instead of 2 and the interior edge have 0 constraints instead of 2, so that the solution global vector has a size of 8?

Shouldn't the layout be something like the following?

(   0): dim 2 offset 0
(   1): dim 2 offset 0
(   2) dim  2 offset   4 constrained 0 1
(   3) dim  2 offset   6 constrained 0 1
(   4) dim  2 offset   8 constrained 0 1
(   5) dim  2 offset  10 constrained 0 1
(   6) dim  4 offset  12 constrained 0 1 2 3
(   7) dim  4 offset  16
(   8) dim  4 offset  20 constrained 0 1 2 3
(   9) dim  4 offset  24 constrained 0 1 2 3
(  10) dim  4 offset  28 constrained 0 1 2 3

Note: P2 with this mesh ends up with all DOF constrained, but I think the middle edge should not be constrained.

Current residual

(0)  0.  0.
(1)  0.  0.
(2)  0.  0.
(3)  0.  0.
(4)  0.  0.
(5)  0.  0.
(6)  0.  0.  0.45  0.6375
(7)  0.  0.  0.  0.
(8)  0.  0.  0.4125  0.45
(9)  0.  0.  -0.4125  -0.45
(10) 0.  0.  -0.45  -0.6375

[knepley]

Yes, to all questions. So we need to look at the code that implements constraints. Where is it in PyLith?

[baagaard-usgs]

MMSTest::testJacboianTaylorSeries() failing due to negative mesh point.

This appeared after @knepley recent push to PETSc knepley/pylith.

cd mmstests/elasticity
make check

[baagaard-usgs]

Negative mesh point was due to uninitialized variable in dmplexsnes.c at line 2182.

P3 test is still failing.

Without Matt's P3 fix, in pylith::fekernels::IsotropicLinearElasticityPlaneStrain::g1v the solution that is passed in is correct for the coordinates.

With Matt's P3 fix, the solution is NOT correct for the coordinates.

Solution for uniform strain

>>> def disp(x, y):
    exx = 0.1
    eyy = 0.25
    exy = 0.3
    disp_x = exx*x + exy*y
    disp_y = exy*x + eyy*y
    return (disp_x, disp_y)

Solution and coordinates in kernel for case when solution matches coordinates:

Thread 1 "test_isotropicl" hit Breakpoint 2, pylith::fekernels::IsotropicLinearElasticityPlaneStrain::g1v (dim=2, numS=1, numA=3, sOff=0x55555617a6c0, sOff_x=0x55555617ad20, s=0x555556187bb0, s_t=0x0, 
    s_x=0x55555618d830, aOff=0x555555dca9f0, aOff_x=0x555555dbc190, a=0x555555b73aa0, a_t=0x0, a_x=0x555555b76360, t=0, x=0x55555618dea0, numConstants=0, constants=0x0, g1=0x55555618e500)
    at /home/brad/src/cig/pylith/libsrc/pylith/fekernels/IsotropicLinearElasticity.cc:61
61      assert(_dim == dim);
(gdb) p x[0]@2
$1 = {-0.86906601088997104, -0.88579160777096477}
(gdb) p s[0]@2
$2 = {-0.35264408342028686, -0.48216770520973284}

[baagaard-usgs]

Created new issue (https://github.com/geodynamics/pylith/issues/121) for higher order discretization order.