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

Thu Feb 28 11:58:39 PST 2008

On Feb 25, 2008, at 11:41 PM, Brad Aagaard wrote:

> 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

I also think this seems like an error, either in the mesh or in the  
solution.  The solution should definitely be converging to the  
analytical solution.

>
>
>> (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).

Yes.  I believe that you are generally better off using quad4/hex8  
elements, as you get better accuracy for the amount of computation  
required.  This is certainly true for linear elements, but we need to  
see what happens with higher-order elements.

>
>
>> (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.

It might be useful to superimpose contours of the computed strain  
field on something that shows the mesh.  Unfortunately, unless you are  
pulling the developers' version, you probably can't output strains  
right now.

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

I haven't tried this, but using something like ParaView, I'm pretty  
sure you could use the calculator to add the computed displacements to  
the coordinates, and then plot the distorted mesh using those.

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

This is definitely something we want in the near future.  There are  
actually two separate possibilities:
1.  A complete restart file, which would contain all information.
2.  The ability to specify initial stresses, which would allow us to  
use either computed or specified stresses to start a simulation.

Thanks for all your work looking at code performance, etc.  This will  
actually be very useful and I'm sure we could use some of this work in  
our next workshop.  Also, at some point I may have some time to play  
around with your meshes to see what I come up with.  If you can post  
them somewhere (e.g., on the CIG website), I'll see what I can come up  
with.

Thanks,
Charles

>
>
> Brad
>

Charles A. Williams
Dept. of Earth & Environmental Sciences
Science Center, 2C01B
Rensselaer Polytechnic Institute
Troy, NY  12180
Phone:    (518) 276-3369
FAX:        (518) 276-2012
e-mail:    willic3 at rpi.edu

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://geodynamics.org/pipermail/cig-short/attachments/20080228/f8ee0115/attachment.htm 