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

[ZHOU Yu]

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://geodynamics.org/pipermail/cig-short/attachments/20130718/5a9d7872/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: attached.tar.gz
Type: application/octet-stream
Size: 2393921 bytes
Desc: not available
URL: <http://geodynamics.org/pipermail/cig-short/attachments/20130718/5a9d7872/attachment-0001.obj>