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

[Brad Aagaard]

On 07/23/2013 08:34 AM, Yu ZHOU wrote:
> Dear sir:
>
> Thank you very much for your instantaneous replies last Friday. I
> redesigned my calculation with considerations of your suggestions and
> did lots of numerical experiments in last few days. I think I've
> figured out this problem.
>
> First, I constructed a purely linear slip distribution that tapers
> into zero at 3 buried edges as Brad suggested, so that FEM nodes can
> exactly simulate the distribution. For Okada solution, the 40kmx12km
> fault is divided into 1kmx1km patches, which I think is an adequate
> resolution but the time consumption is about 1 hour. The slip on a
> particular patch is the average value of its four vertices. The
> discrepancy is cut down greatly for the slip distribution is
> "natural" for FEM, but it still exists.

For the Okada solution you can use variable patch sizes to minimize the 
total number of patches. For example, in the region where the slip is 
uniform, you only need 1 patch. In most cases this type of optimization 
results in using far fewer patches compared with using a uniform patch size.

> The discrepancy comes from 3 factors in addition to "unnatural" slip
> distributions: 1) zero disp boundary conditions 2) element size near
> fault, Sf, and that near exterior boundary, Se 3) linear nature of
> shape functions of tetrahedron elements
>
> 1.Matthew suggested using Okada solution results on the boundary as
> consistent Dirichlet conditions, but I think it's only helpful for
> developers, not users, because people choose FEM to simulate
> coseismic  deformation in heterogeneous lithosphere or with
> topography, not in homogeneous elastic half space. In almost all
> published papers about coseismic and postseismic simulation with FEM,
> boundary conditions are set to zero displacement. In fact, however,
> displacement can never be absolutely zero, so zero disp bd condition
> is inconsistant. As a result, the closer to the exterior boundary,
> the larger the discrepancy is, so FEM results are sufficiently
> accurate only within a domain far from boundary and close to the
> fault rupture. In my example, for a 2400kmx2400kmx1200km mesh volume,
> the discrepancy is below 2% in the domain where |x|<300km and
> |y|<300km. A larger accurate domain requires larger mesh volume
> dimensions. I suggest this fact should be mentioned in Pylith
> manual.

The target accuracy in most applications is to keep the numerical errors 
well below the uncertainty in the observations. Because observations 
usually have uncertainties that are independent of the deformation (for 
example uncertainties in GPS displacements are often in the mm range), 
for practical problems the absolute error is more important than 
relative error (which you seem to be using as an error metric). We have 
not had difficulty in setting up problems to meet this target accuracy. 
Remember that the earth is not an infinite medium and slip is not 
discontinuous, so the finite-element solution may involve fewer 
approximations than the Okada solution.

> 2.Usually element size near the fault is very fine, gradually coarser
> to boundaries. Certainly smaller Se enlarges the accurate domain, but
> element number will soar.  It's very interesting that accurate domain
> doesn't always expand with smaller Sf. In my example, accurate domain
> is largest when Sf is 3km if it loops over [1km, 2km, 3km, 4km].

You may not have an optimal spatial variation in the discretization 
size. P1 basis functions can capture a linear variation in the 
displacement field, so a good rule of thumb is that the discretization 
size should be proportional to the gradient in the strain. This means 
you need the finest discretization size in the regions where there is a 
change in the gradient of the slip or displacement field.

> 3. Shape function of hex element may be nonlinear in Pylith(I cannot
> find expressions of shape functions in Pylith manual), but it's not
> convenient to generate variable sized hex mesh in CUBIT in
> complicated geometry compared with tet, and I have not do any tests
> with hex meshes. Perhaps hex meshes behave much better with the same
> number of nodes?

You can see exactly what basis functions are being used by turning on 
the FIATSimplex journal as in examples/3d/tet4/pylithapp.cfg (line 16). 
Setting the journal flag to True or 1 will display the basis functions 
and their derivatives evaluated at the quadrature points in addition to 
the quadrature information. Hex cells with P1 basis functions do contain 
more terms than tet cells with p1 basis functions. As a result, for a 
given discretization size, the hex cells tend to perform better. As you 
note, it is more difficult to vary the discretization size with hex8 
cells, so which cell to use depends on the problem.

Regards,
Brad