<div dir="ltr">On Thu, Jul 18, 2013 at 5:52 AM, ZHOU Yu <span dir="ltr"><<a href="mailto:zhouyuking@foxmail.com" target="_blank">zhouyuking@foxmail.com</a>></span> wrote:<br><div class="gmail_extra"><div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Dear Pylith developers:<br><br>I build up a homogeneous FEM model to compare coseismic ground displacement field output by Pylith with Okada analytical solution, <br>
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.<br>I would be very grateful if you can help me with this discouraging problem.<br>
<br>My FEM model settings are as the following:<br><div style="text-align:left;margin-left:80px"><font size="1"><font><span style="font-family:Times New Roman">Model dimension: 1800kmx2000kmx600km</span><br><span style="font-family:Times New Roman">Fault geometry: strike=0, dip=90, 40km long, 12km wide</span><br>
<span style="font-family:Times New Roman">Tetrahedron element size: 1km on or near the fault, gradually larger to 100km on the boundaries.</span><br><span style="font-family:Times New Roman">Coseismic fault slip: 10m strike slip, homogeneous on the fault</span><br>
<span style="font-family:Times New Roman">Material type: homogeneous Poisson solid</span><br><span style="font-family:Times New Roman">Boundary conditions: fixed along the normal direction for 4 lateral surfaces and bottom surface, free on ground surface</span></font><br>
</font></div><br>In the attached tarball please find <span style="font-family:Times New Roman">mesh_shortfaullt.jou</span>, which is the journal file for mesh generating in CUBIT. Its output data file <span style="font-family:Times New Roman">mesh_shortfault.exo</span> is not included due <br>
to its big size(>50Mbyte). All other files that Pylith needs can also be found in the tarball, please just run:<br><div style="margin-left:80px"> <span style="font-family:Times New Roman">pylith eqsim.cfg</span><br></div>
Output hdf5 files will be stored in <span style="font-family:Times New Roman">./output</span>. Then Please run the matlab script <span style="font-family:Times New Roman">homoslip_diff.m</span>. <br>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.<br>
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<br> 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.<br>
</blockquote><div><br></div><div>This is a fundamental point about the meaning of FEM solutions. In order to compare to another function, you want</div><div><br></div><div>  error = \int_\Omega |u - u^pylith|^2</div><div>
    = \sum_ele \int_ele |u - u^pylith|^2</div><div>    = sum_ele \sum_q |u(x_q) - \sum_basis u^pylith_b \psi_b(x_q)|^2 w_q</div><div><br></div><div>which is just the evaluation of a functional over the FE space. This error should converge to zero like</div>
<div>as h^{-1}, or have slope 1 on a log graph.</div><div><br></div><div>Second, have you input the Okada solution and looked at the residual? This is a good way to tell if the BC you impose</div><div>are consistent with your exact solution.</div>
<div><br></div><div>    Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">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, <br>
the result didn't improve as much.<br><br>Then I thought about interpolation method. Since variables in a tetrahedron are linear, it should be right to employ bilinear interpolation.<br><br>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<br>
 near the fault(see line 41 in <span style="font-family:Times New Roman">mesh_shortfaullt.jou</span>). The number of nodes jumped to bigger than 500,000, but such an excessively dense mesh didn't help either.<br><br>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<br>
dense, and slip distribution is uniform on the entire fault surface, this effect will be neglectable.<br><br>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.<br>
<br><br><br><br>_______________________________________________<br>
CIG-SHORT mailing list<br>
<a href="mailto:CIG-SHORT@geodynamics.org">CIG-SHORT@geodynamics.org</a><br>
<a href="http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short" target="_blank">http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short</a><br></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>