[CIG-SHORT] Meshing Schemes, Mixed Meshes, Mesh Output, & Initial Conditions

[Brad Aagaard]

Ravi Kanda wrote:
> 1) MESHING RESOLUTION TESTS (MeshAccuracy1.pdf):
> (a) Quad4 meshes: The "Best" model (model 1) I can come up with has a 
> fixed interval (size, ds = 20 km).  For Quad4 elements, overall error 
> decreases steadily from ds = 80 km to ds = 20 km (as mesh resolution 
> gets finer), but when it reaches ds = 10 km, the error is huge.  Is this 
> due to round-off errors - i.e. too many mesh elements for the applied 
> fault displacement (10 cm) and the resulting strain field?  I am not 
> sure I understand what is going on here...

It looks like there must be an error in the ds = 10km run. The solution 
should converge. Did the solution converge in the PETSc solve? Did it 
hit the maximum number of iterations that you specified? As discussed in 
the PyLith manual, use the following PETSc options (--petsc.XXX) to 
monitor the solution convergence.

ksp_monitor = true
ksp_view = true

> (b) Tri3 meshes: What I am finding is that overall, a triangular mesh 
> gives me poorer results than a quadrilateral mesh having the same CUBIT 
> mesh interval size (ds = 20 km, model 1).  I realize there is no formal 
> mathematical comparison between unstructured Quad4 and Tri3 meshes. 
 > ...   I am not sure why a triangular mesh
> is giving me poorer performance in general, given it fits the geometry 
> much better.

There is a significant difference in the accuracy of tri3 versus quad4 
and tet4 versus hex8 cells due to the terms in the shape functions. Tri3 
and tet4 cells can exactly represent a linear variation in the 
displacement field, so they have constant strain within the cell. Quad4 
and hex8 cells have more terms in the basis functions so it can better 
capture more complex variations. This is discussed in most 
finite-element texts.

We showed this on the AGU PyLith poster (I will post it on the CIG 
website on Tue) for the strike-slip benchmark. The integrated error with 
1 km resolution hex8 cells was nearly the same as with 500m resolution 
tet4 cells. We obtained nice convergence with each type of cell (the 
error was always less with the hex8 cells at the same resolution).

> (c) Adaptive meshes: Adaptive meshs also give poorer results compared to 
> model 1.  I illustrate this with two different meshes having variable 
> resolution. Again, I realize they cannot me mathematically compared in a 
> formal sense to a "regular interval" mesh, but I don't understand why 
> the adaptive mesh should give a larger error, given almost all the 
> strain is concentrated, where I have as good a resolution (on average) 
> as in model 1.
> 
> So, I want to know whether Pylith has been tested with adaptive meshes, 
> and what kind of mesh resolution tests it has been subjected to.

Remember it is not the magnitude of the displacements that determine 
what resolution you need but the gradient in strain. A uniform strain 
field can be modeling with a single cell, but a large gradient in strain 
requires very small cells. The gradient in strain may be decaying much 
slower with distance from the fault than the displacements, so you may 
be coarsening too fast with distance.

We have not done testing with adaptive meshes. Designing meshes for a 
particular problem is not a software development issue but an 
application issue. Generaly, only when one wants to do AMR does it enter 
into the software development. I have been surprised by how little 
interest there has been within our working group for generating optimal 
meshes. Only Charles Williams, Charles Norton, and I from the short-term 
tectonics community attended the AMR workshop and no one responded to 
the call for participation in using deal.ii to experiment with AMR for 
tectonics problems.

> 2) MIXED-MESHES (MeshAccuracy2.pdf):
> Based on the Tri3 vs. Quad4 element performance above, I created a model 
> geometry where the wedge-shaped region is a different volume - and there 
> are four materials in this model.  I am having a hard time making CUBIT 
> match the mesh across two of the volumes ("purple" & "red" in attached 
> PDF) - the nodes simply don't overlap.  I tried both TRI3 & QUAD4 
> elements.  Would the resulting non-coincident nodes at the material 
> interfaces be a problem for Pylith?  I am still trying to figure out 
> what combination of interval matching (size vs. number) can overcome this.

Mixing cells types is not supported and we don't have any plans to do 
so. Non-coincident nodes at the material interfaces is a problem. This 
is a "bad" mesh. I am glad PyLith spit out an error message, but it 
should have been caught on importing the mesh.

> 3) OUTPUT DISTORTED MESH FROM PYLITH RUN:
> Another issue I am facing is that even thought my problem is fully 
> elastic, applying significant slip (say on the order of several 100 m 
> to  1 km, which is still much less than smallest element size of 20 km 
> in model 1) on the fault to simulate long term motion would result in 
> the distortion of the fault surfaces as well as the free surface at the 
> top.  One way to capture this would be to have access to the distorted 
> mesh at the end of the pylith run, along with the displacements VTK 
> file.  This could then be used as an input for the next stage.
> At this stage, is it possible to hack Pylith to export the distorted 
> mesh?  If not, are there any plans for Pylith to output the mesh in the 
> near future?

I am not sure that I follow what you are asking. PyLith currently uses 
an infinitesimal strain formulation, so the "distorted" mesh is simply 
the original mesh with the nodes displaced by the solution field. The 
original mesh is obviously available as is the solution.

> 4) INITIAL CONDITION FOR PYLITH RUN:
> In order to accomplish (3) above, in addition to the distorted mesh, I 
> also need to be able to set the initial condition to the previous run's 
> output displacement.  So, again, is there some way to make Pylith accept 
> a prescribed initial condition?  If not, is this something that is in 
> the plans for the near future?

You can update the position of the vertices in the mesh with the 
solution to create a new mesh and create a spatial database file with 
whatever initial condition you want (e.g., computed from the end of a 
previous solution, etc). Once we get HDF5 output going (currently 
scheduled for version 1.2), this could be done with only a few lines of 
Python code.

Brad