[CIG-SHORT] significant discrepancy between coseismic displacement field of homogeneous FEM model and Okada analytical solution

[Matthew Knepley]

On Thu, Jul 18, 2013 at 4:42 PM, Brad Aagaard <baagaard at usgs.gov> wrote:

> Yu Zhou,
>
> The main issue in comparing PyLith prescribed slip results with Okada is
> that Okada uses uniform slip patches while PyLith uses a continuous slip
> field that is C0 continuous. At the edges of a large uniform slip patch,
> there is a stress concentration, which would require very fine
> discretization with PyLith. If you want to compare PyLith against an Okada
> solution, the preferred approach is to construct a slip distribution that
> tapers linearly to zero at the buried edges of the fault. This is more
> difficult to represent in Okada. We have used many overlapping uniform slip
> patches of varying sizes stacked together to form a slip distribution with
> very fine stair stepping to approximate the linear variation represented
> exactly in PyLith. Our results show nice convergence using the error metric
> Matt mentioned.
>

This is probably worth mentioning in the manual since its is the first
thing many people are going to try. Better yet would be a benchmark
paper with the other codes that do viscoelasticity and faults. Then we
could just point people to the paper with everything worked out.

    Matt


> Regards,
> Brad
>
>
>
> On 07/18/2013 03:52 AM, ZHOU Yu wrote:
>
>> Dear Pylith developers:
>>
>> I build up a homogeneous FEM model to compare coseismic ground
>> displacement field output by Pylith with Okada analytical solution,
>> but I find unneglectable discrepancy between the two.  I've examined
>> possible reasons and tried hard, but I cannot eliminate it through all my
>> efforts.
>> I would be very grateful if you can help me with this discouraging
>> problem.
>>
>> My FEM model settings are as the following:
>> Model dimension: 1800kmx2000kmx600km
>> Fault geometry: strike=0, dip=90, 40km long, 12km wide
>> Tetrahedron element size: 1km on or near the fault, gradually larger to
>> 100km on the boundaries.
>> Coseismic fault slip: 10m strike slip, homogeneous on the fault
>> Material type: homogeneous Poisson solid
>> Boundary conditions: fixed along the normal direction for 4 lateral
>> surfaces and bottom surface, free on ground surface
>>
>>
>> In the attached tarball please find mesh_shortfaullt.jou, which is the
>> journal file for mesh generating in CUBIT. Its output data file
>> mesh_shortfault.exo is not included due
>> to its big size(>50Mbyte). All other files that Pylith needs can also be
>> found in the tarball, please just run:
>>   pylith eqsim.cfg
>>
>> Output hdf5 files will be stored in ./output. Then Please run the matlab
>> script homoslip_diff.m.
>> This Matlab script reads displacement at every node on the ground surface
>> in FEM mesh, then interpolates disp field at a group of gridded points.
>> Theoretical disp field is also calculated with Okada solution at those
>> gridded points. Differences between Okada disp and FEM disp are shown as a
>> vector map
>>   and 3 grid images. Discrepancy for disp-Y component, which is along
>> strike direction, is larger than 6% in most area, increasing with distance
>> from the fault.
>>
>> At the beginning, the interval size was set to 2km on fault interface, so
>> I thought perhaps the mesh was not dense enough near the fault. After I
>> changed it to 1km,
>> the result didn't improve as much.
>>
>> Then I thought about interpolation method. Since variables in a
>> tetrahedron are linear, it should be right to employ bilinear interpolation.
>>
>> Later, I thought that strain component exy is very large due to pure
>> strike slip in the neighbourhood of the fault, I set interval size to 0.5km
>> along x direction
>>   near the fault(see line 41 in mesh_shortfaullt.jou). The number of
>> nodes jumped to bigger than 500,000, but such an excessively dense mesh
>> didn't help either.
>>
>> I've also considered that slip at a single node in FEM means linear slip
>> distribution tapering into zero at the neighbouring nodes. However, if the
>> mesh is sufficiently
>> dense, and slip distribution is uniform on the entire fault surface, this
>> effect will be neglectable.
>>
>> After all these computation and consideration, now I wonder whether this
>> discrepancy is kind of intrinsic error of Pylith software and cannot be
>> eliminated by users.
>>
>>
>>
>> ______________________________**_________________
>> CIG-SHORT mailing list
>> CIG-SHORT at geodynamics.org
>> http://geodynamics.org/cgi-**bin/mailman/listinfo/cig-short<http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>>
>>
> ______________________________**_________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://geodynamics.org/cgi-**bin/mailman/listinfo/cig-short<http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://geodynamics.org/pipermail/cig-short/attachments/20130718/8325ed4f/attachment-0001.html>