# How to impose Neumann shear traction on a diagonal boundary? **URL:** \ **Category:** PyLith\ **Created:** [April 1, 2020, 8:13pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187 "2020-04-01T20:13:53Z")\ **Posts on this page:** 13\ **Page:** 1
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 1, 2020, 8:13pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/1 "2020-04-01T20:13:53Z")
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 ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [April 1, 2020, 9:25pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/2 "2020-04-01T21:25:50Z")
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). ---
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 2, 2020, 9:10pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/3 "2020-04-02T21:10:11Z")
![c_view](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/1X/6ed26880218839e55afedb6adcf942781c6dc840.png) 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](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](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} } ---
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 7, 2020, 7:02pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/4 "2020-04-07T19:02:32Z")
[.json file](https://drive.google.com/file/d/1ZhfpbJXRZr_8v0Ju2bhKWmU1aIJyExGC/view?usp=sharing) Here is a link to the .json file on my google drive. You should be able to view it. ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [April 7, 2020, 7:41pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/5 "2020-04-07T19:41:13Z")
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. ```auto output.cell_info_fields = [initial_value] ``` ---
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 8, 2020, 6:46pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/6 "2020-04-08T18:46:11Z")
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 ---
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 8, 2020, 8:43pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/7 "2020-04-08T20:43:50Z")
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 ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [April 9, 2020, 1:59am UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/8 "2020-04-09T01:59:10Z")
What solver (PETSc) parameters are you using? What does output (please provide the entire log) show when you use: ```auto [petsc] ksp_monitor = true ksp_converged_reason = true ``` ---
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 9, 2020, 6:42pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/9 "2020-04-09T18:42:50Z")
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](https://drive.google.com/file/d/16dpr1mXQmgN2nXb5BwQcLSVjvROdidO_/view?usp=sharing) Joe ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [April 9, 2020, 9:56pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/10 "2020-04-09T21:56:05Z")
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: ```auto [pylithapp.petsc] pc_type = lu ``` For production runs, we recommend using the algebraic multigrid preconditioner: ```auto 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`. ---
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 13, 2020, 8:25pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/11 "2020-04-13T20:25:48Z")
![model_c2](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/1X/b09da8f71b7f4d09ed0b7c45eacf41f92e6b74af.png) 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 ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [April 13, 2020, 8:43pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/12 "2020-04-13T20:43:12Z")
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? ---
**Author:** ![RipUMN](https://avatars.discourse-cdn.com/v4/letter/r/57b2e6/32.png) [@RipUMN](https://community.geodynamics.org/u/RipUMN)\ **Post date:** [April 13, 2020, 8:49pm UTC](https://community.geodynamics.org/t/how-to-impose-neumann-shear-traction-on-a-diagonal-boundary/1187/13 "2020-04-13T20:49:41Z")
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.