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 size="2"><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><br>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>