How to impose Neumann shear traction on a diagonal boundary?

PyLith
1

I’m trying to impose shear on a diagonal boundary of my model domain (approximating a subduction zone). I have successfully run a simple model with the same Neumann shear BC on a vertically oriented boundary, but I can’t seem to get it to work on the diagonal, megathrust, boundary. I have tried:
– Neumann BC on the diagonal boundary
– Tilting the model domain such that the diagonal boundary is horizontal, then altering the gravity orientation appropriately, then imposing the Neumann BC
– Simplifying the mesh

In all cases the model will not converge (I have tried easing the convergence parameters too). Please help. Can I approach this problem differently? Am I missing something? Please let me know.

Joe

2

Applying a Neumann BC to a diagonal boundary should work. Most likely there is something else wrong with the problem setup. It is always helpful if you supply a diagram/cartoon illustrating the problem you are trying to solve. Please also include the PyLith version information (we prefer the output of pylith --version). If you need help debugging your .cfg files, then include the .json parameter file from the simulation in question (preferably one that matches your diagram).

1 Like
3

c_view
c_view821×632 13.9 KB

Hi! Thank you for replying. I have attached a graphic of the model domain I’m using. It’s basically a rectangular solid halved by a diagonal that includes 3 material types. I want to impose a shear traction along the diagonal side to approximate the shear stress of a subducting plate. Ideally I want to vary the stress by depth, but right now I’m just trying to get it to run. Also, if it matters, once I get it to run I’m going to run the models on an 8 processor cluster, not my laptop.

I’m not allowed to attach documents to this reply because I’m a new member, so I’m going to paste below the text output of pylith --version, and my .json parameter file. **update on this: the .json file was too large to paste into this reply. I would like to provide it to you. Can you update my user status so I can attach it, or provide another method for me to send you the file? **

Thanks again! I really do appreciate any help you can offer.

Joe

Pylith version (output of pylith --version)

Platform:
Hostname: Josephs-MacBook-Air-2.local
Operating system: Darwin
Kernel: 19.3.0
Version: Darwin Kernel Version 19.3.0: Thu Jan 9 20:58:23 PST 2020; root:xnu-6153.81.5~1/RELEASE_X86_64
Machine: x86_64
Processor: i386

PyLith
Release v2.2.1

PETSc
Configured on 2017-08-03 14:01:44 -0500, GIT branch: branch-not-available, revision: v3.7.6-4826-gd686aaf

Spatialdata
Release v1.9.10

MPI standard: 3.0, implementation: MPICH, version: 3.1.3
HDF5 version: 1.8.11
NetCDF version: 4.4.1.1
Proj version: 480

Python
v2.7.16 of CPython compiled with GCC 4.2.1 Compatible Apple LLVM 11.0.0 (clang-1100.0.32.4) (-macos10.15-objc-s
FIAT: v0.9.9 from /Users/JoeRippke/Pylith/lib/python2.7/site-packages/FIAT
spatialdata: v1.9.10 from /Users/JoeRippke/Pylith/lib/python2.7/site-packages/spatialdata
netCDF4: vnot found from –
h5py: v2.2.1 from /Users/JoeRippke/Pylith/lib/python2.7/site-packages/h5py
pyre: v0.8.1 from /Users/JoeRippke/Pylith/lib/python2.7/site-packages/pythia-0.8.1.18-py2.7.egg/pyre
numpy: v1.11.2 from /Users/JoeRippke/Pylith/lib/python2.7/site-packages/numpy-1.11.2-py2.7-macosx-10.11-intel.egg/numpy

If you publish results based on computations with PyLith please cite the following:
(use --include-citations during your simulations to display a list specific to your computation):

@Manual{PyLith:software,
title = {PyLith v2.2.1},
author = {Aagaard, B. and Knepley, M. and Williams, C.},
organization = {Computational Infrastructure for Geodynamics (CIG)},
address = {University of California, Davis},
year = {2017},
doi = {http://doi.org/10.5281/zenodo.886600}
}

@Manual{PyLith:manual,
title = {PyLith User Manual, Version 2.2.1},
author = {Aagaard, B. and Knepley, M. and Williams, C.},
organization = {Computational Infrastructure for Geodynamics (CIG)},
address = {University of California, Davis},
year = {2017},
note = {http://www.geodynamics.org/cig/software/pylith/pylith_manual-2.2.1.pdf}
}

@Article{Aagaard:Knepley:Williams:JGR:2013,
author = {Aagaard, B.~T. and Knepley, M.~G. and Wiliams, C.~A.},
title = {A domain decomposition approach to implementing fault slip in finite-element models of quasi-static and dynamic crustal deformation},
journal = {Journal of Geophysical Research Solid Earth},
year = {2013},
volume = {118},
pages = {3059–3079},
doi = {10.1002/jgrb.50217}
}

4

.json file

Here is a link to the .json file on my google drive. You should be able to view it.

5

I don’t see any obvious error in your .cfg file. Some things to check:

  1. Do you have any duplicate surfaces in your geometry (so that the materials are actually disconnected)? You should use imprint all and merge all at the end of your geometry creation in CUBIT/Trelis.

  2. Check your node sets in CUBIT/Trelis to make sure they match what you think they should be.

  3. Output the initial values of the fault tractions to make sure your spatial database is providing the values that you think it is. See examples/3d/hex8/step02 for an example. Note that the components will be tangential (horizontal), tangential (other = normal x tangential horiz), and normal.

This turns on output of the initial values within the Neumann BC settings.

output.cell_info_fields = [initial_value]
6

Thank you! I’m glad to hear there’s nothing obviously wrong with my .cfg files. Although, it would be nice if there was an easy fix. I’ll check some things out as you suggest.

Joe

7

I verified that my mesh has no duplicate surfaces and that my node sets are correct, and I reviewed my initial traction values as you suggested. I improved my initial traction, so I’m happy about that, but the model still won’t converge. I’m going to keep at it. Maybe I’ll figure it out. Please let me know if anything else occurs to you.
Thanks again,
Joe

8

What solver (PETSc) parameters are you using? What does output (please provide the entire log) show when you use:

[petsc]

ksp_monitor = true
ksp_converged_reason = true
9

PETSc parameters:

[pylithapp.petsc]
malloc_dump =

ksp_rtol = 1.0e-8
ksp_atol = 1.0e-12
ksp_max_it = 200
ksp_gmres_restart = 50

ksp_monitor = true
ksp_converged_reason = true
ksp_error_if_not_converged = true
log_view = true

Here’s the entire output log:
output log

Joe

10

It looks like the solver is getting stuck (the residual isn’t decreasing). This could be solver settings or you have some very distorted cells (check your mesh quality).

It looks like you are using the default PETSc linear solver and preconditioner settings. In addition to specifying the tolerances, we recommend setting the solver and/or preconditioner.

No fault

For debugging purposes, we recommend the LU preconditioner:

[pylithapp.petsc]
pc_type = lu

For production runs, we recommend using the algebraic multigrid preconditioner:

pc_type = ml

See the PyLith manual and online tutorials for more information on solver settings. When you have a fault, we recommend the custom fault preconditioner + algebraic multigrid preconditioning. These settings are in share/settings/solver_fault_fieldsplit.cfg.

11

model_c2
model_c21075×683 22.7 KB

I found a solution!

I blunted the pointy tips of my model domain (see the attached image, and compare with the image in a previous post). Now the model runs and produces output that makes sense.

I’m not sure why this is the fix.

I combed through my code so many times, and I adjusted my solver settings many times, but nothing worked until I changed the shape of my model domain.

Thanks again for your help!
Joe

12

My guess is the issue is related to an overlap in the Dirichlet and Neumann BC at the corners. With this change in the domain, do your Dirichlet and Neumann BC now have zero points in common whereas before they did?

13

Should be about the same. I imposed no slip dirichlet BCs on the blunted tips of this new domain. So, the neumann condition is still bordered by dirichlet BCs on all sides. Although, my original domain had a free surface dirichlet BC on the top surface which bordered the neumann BC at a point. That being said, blunting the point at the surface didn’t fix the problem, I had to blunt both points for it to work.