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

[Matthew Knepley]

On Thu, Jul 18, 2013 at 5:52 AM, ZHOU Yu <zhouyuking at foxmail.com> 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.
>

This is a fundamental point about the meaning of FEM solutions. In order to
compare to another function, you want

  error = \int_\Omega |u - u^pylith|^2
    = \sum_ele \int_ele |u - u^pylith|^2
    = sum_ele \sum_q |u(x_q) - \sum_basis u^pylith_b \psi_b(x_q)|^2 w_q

which is just the evaluation of a functional over the FE space. This error
should converge to zero like
as h^{-1}, or have slope 1 on a log graph.

Second, have you input the Okada solution and looked at the residual? This
is a good way to tell if the BC you impose
are consistent with your exact solution.

    Matt


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



-- 
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/41cc4d42/attachment.html>