[CIG-SHORT] Spherical 3d gravity initial stress
Brad Aagaard
baagaard at usgs.gov
Wed Nov 9 09:04:45 PST 2016
Demian,
Your header shows the grid is
num-x = 33
num-y = 21
num-z = 28
so the total number of points should be 33*21*28, but I only see 592
points in the file.
Brad
On 11/09/2016 08:55 AM, Demian Gomez wrote:
> Brad,
>
> I'm analyzing the far-field (> 800 km from a thrust fault) co-seismic
> displacements and the effect of the asthenospheric stiffness in
> the lithospheric stress guide. It's mostly a "theory" analysis and I
> cannot neglect the effect of sphericity since this would be highly
> criticized by reviewers.
>
> I am trying to implement a spatial grid db as you suggested, but I am
> having a few problems. First for all, for some reason Pylith is crashing
> with my spatial grid with the following error:
>
> RuntimeError: Error occurred while reading spatial database file
> 'spatialdb/initial_stress-layer1.spatialdb'
> I/O error while reading SimpleGridDB data.
>
> I thought the problem could be the size of the spatialdb so I reduced
> the number of
> points but the problem persisted.
>
> I then downloaded a traction db example from a benchmark
> (https://github.com/geodynamics/pylith_benchmarks/tree/master/dynamic/scecdynrup/tpv102
> <https://github.com/geodynamics/pylith_benchmarks/tree/master/dynamic/scecdynrup/tpv102>)
> to compare it against my file. I couldn't see any remarkable differences
> (except the coordinate system, number of points, etc) so I used my
> header (modifying num-x, num-y and num-z) with the data portion of the
> benchmark spatial db. Even though the number of values (columns) in this
> spatial db do not agree with the 6 stress components that I need (the
> benchmark uses 7 values) and the values itself are completely
> meaningless in terms of stress, Pylith does now pass the material
> initialization. It then crashes with the second layer where I left the
> spatial db unmodified.
>
> There seems to be some kind of problem with my data portion of the
> spatial db, but I can't find it. I looked for CRLF/LF problems, field
> sizes, etc just in case, but found no differences. Attached you will
> find the spatial grid db that crashes Pylith. I'd appreciate if you
> could take a look at it to see if there is anything obviously wrong.
>
> Demián
>
> --
> *Dr. Demián D. Gómez*
> Postdoctoral Researcher
> The Ohio State University - School of Earth Sciences
> 275 Mendenhall Laboratory
> 125 South Oval Mall
> Columbus, Ohio 43210
> Cell: +1 (901) 900-7324 <tel:%2B1%20%28901%29%20900-7324>
> email: gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>
>
>
> On Wed, Nov 9, 2016 at 9:50 AM, Brad Aagaard <baagaard at usgs.gov
> <mailto:baagaard at usgs.gov>> wrote:
>
> Demian,
>
> The code for initializing the initial stress and strain fields is in
> the method ElasticMaterial::_initializeInitialStress() in
> libsrc/pylith/materials/ElasticMaterial.cc.
>
> Before you head in that direction, I would seriously consider the
> size of the approximation in using a flat domain compared to a
> spherical domain relative to the other assumptions in your problem.
> In other words, do your have sufficient resolution in all of your
> other parameters that you know you need to use a domain with a
> spherical cap?
>
> Regards,
> Brad
>
>
>
> On 11/08/2016 04:08 PM, Demian Gomez wrote:
>
> Hi Brad,
>
> Thanks for your answer.
>
> Would it be possible to modify Pylith's code to calculate the
> initial
> tractions during runtime without having to specify them through a
> spatial database? I can try to do it if you can provide me with some
> orientation on where I should start looking. It seems that for a
> model
> constructed in ECEF coordinates, calculating stress_xx,
> stress_yy and
> stress_zz should be trivial if I have access to the cell
> centroids and
> the material density.
>
> Thanks again,
> Demián
>
> --
> *Dr. Demián D. Gómez*
> Postdoctoral Researcher
> The Ohio State University - School of Earth Sciences
> 275 Mendenhall Laboratory
> 125 South Oval Mall
> Columbus, Ohio 43210
> Cell: +1 (901) 900-7324 <tel:%2B1%20%28901%29%20900-7324>
> <tel:%2B1%20%28901%29%20900-7324>
> email: gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>
> <mailto:gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>>
>
>
>
> On Mon, Nov 7, 2016 at 12:42 PM, Brad Aagaard <baagaard at usgs.gov
> <mailto:baagaard at usgs.gov>
> <mailto:baagaard at usgs.gov <mailto:baagaard at usgs.gov>>> wrote:
>
> Demian,
>
> I think you can get things to work with a 1-D spatial
> database with
> some minor changes.
>
> WARNING: This will only work if stress_xx = stress_yy =
> stress_zz
> and stress_xy = stress_yz = stress_xz = 0. This is due to
> the fact
> that the orientation of the gravity vector is different in the
> projected coordinate system than it is in the geocentric
> coordinate
> system.
>
> You should construct your model in a geocentric local Cartesian
> coordinate system. This is a geocentric coordinate system
> with the
> origin at the surface.
>
> The initial tractions should be in a projected (flat) geographic
> coordinate system, so that stress(z) = integral of rho*g
> from 0 to z.
>
> As an aside, the SimpleGridDB provides much faster interpolation
> than SimpleDB for a logically Cartesian grid because it can
> find the
> relevant points without a global search. The points need to
> conform
> to a grid, but the x, y, and z coordinates do not have to be
> spaced
> uniformly.
>
> Regards,
> Brad
>
>
>
> On 11/06/2016 01:16 PM, Demian Gomez wrote:
>
> Dear all,
>
> I am trying to turn on gravity in a quasi-static problem
> that
> uses a 3D
> spherical cap model with layers (see attached figure for
> geometry and
> displacement field without gravity). The coordinate
> system of
> the mesh
> is ECEF and in pylithapp.cfg I've defined coordsys =
> spatialdata.geocoords.CSGeo. My spatial databases for
> material
> properties, etc are defined using:
>
> num-locs = 1
> data-dim = 0
> space-dim = 3
> cs-data = geographic {
> to-meters = 1.0
> space-dim = 3
> ellipsoid = sphere
> datum-vert = ellipsoid
> is-geocentric = true
> } // cs-data
>
> Since the layers are uniform, I use
> db_properties.query_type =
> nearest.
>
> To turn on gravity, I need to set the initial
> lithostatic stress =
> g*rho*h for each element (I can't use a uniform database
> anymore), but
> because this is a large model (~15e6 elements), the
> sizes of the
> initial
> stress databases are > 300 MB which makes Pylith take a very
> long time
> to initialize the materials (so far I've waited almost
> 12 hours and
> Pylith was still initializing the first layer!).
>
> Following this thread:
>
> http://lists.geodynamics.org/pipermail/cig-short/2014-July/001818.html
> <http://lists.geodynamics.org/pipermail/cig-short/2014-July/001818.html>
>
> <http://lists.geodynamics.org/pipermail/cig-short/2014-July/001818.html
> <http://lists.geodynamics.org/pipermail/cig-short/2014-July/001818.html>>
>
> I have now attempted to create an initial stress spatial
> database that
> is only 1D as follows:
>
> #SPATIAL.ascii 1
> SimpleDB {
> num-values = 6
> value-names = stress-xx stress-yy stress-zz stress-xy
> stress-yz stress-xz
> value-units = Pa Pa Pa Pa Pa Pa // units
> num-locs = 2
> data-dim = 1
> space-dim = 3
> cs-data = geo-local-cartesian {
> origin-lon = -72.6425
> origin-lat = -35.6331
> origin-elev = 0
> to-meters = 1.0
> ellipsoid = sphere
> datum-vert = ellipsoid
> } // cs-data
> }
> 0.0000 0.0000 0.0000 0.000
> 0.000
> 0.000 0.00 0.00 0.00
> 0.0000 0.0000 -104911.1753 -2674950661.887
> -2674950661 <tel:2674950661>.887
> -2674950661 <tel:2674950661> <tel:2674950661
> <tel:2674950661>>.887 0.00 0.00 0.00
>
>
> and setting db_initial_stress.query_type = linear.
>
> My idea was that Pylith would map the
> geo-local-cartesian from the
> spatial database to ECEF coordinates in geocoords.CSGeo.
> However, I got
> the following error during execution:
>
> mpinemesis: spatialdb/SimpleDBQuery.cc:352: void
>
> spatialdata::spatialdb::SimpleDBQuery::_findLinePt(std::vector<spatialdata::spatialdb::SimpleDBQuery::WtStruct>*):
> Assertion `locIndexA >= 0' failed.
>
> My questions are:
>
> 1) Is it possible to mix coordinate systems the way I did?
> 2) Is there another way to tell Pylith that the initial
> stresses
> in my
> model vary linearly with depth, even when I have an ECEF
> coordinate
> system? I would like to avoid having to assign the
> initial stress to
> every element one by one (which leads to a very long
> initialization time).
>
> Any suggestions will be greatly appreciated. Let me know
> if you need
> more information.
>
> Cheers,
> Demián
>
> --
> *Dr. Demián D. Gómez*
> Postdoctoral Researcher
> The Ohio State University - School of Earth Sciences
> 275 Mendenhall Laboratory
> 125 South Oval Mall
> Columbus, Ohio 43210
> Cell: +1 (901) 900-7324
> <tel:%2B1%20%28901%29%20900-7324> <tel:%2B1%20%28901%29%20900-7324>
> email: gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>
> <mailto:gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>>
> <mailto:gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>
> <mailto:gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>>>
>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> <mailto:CIG-SHORT at geodynamics.org>
> <mailto:CIG-SHORT at geodynamics.org
> <mailto:CIG-SHORT at geodynamics.org>>
>
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short>>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org <mailto:CIG-SHORT at geodynamics.org>
> <mailto:CIG-SHORT at geodynamics.org
> <mailto:CIG-SHORT at geodynamics.org>>
>
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short>>
>
>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org <mailto:CIG-SHORT at geodynamics.org>
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org <mailto:CIG-SHORT at geodynamics.org>
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
> <http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>
>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>
More information about the CIG-SHORT
mailing list