[CIG-SHORT] Spherical 3d gravity initial stress
Brad Aagaard
baagaard at usgs.gov
Wed Nov 9 06:50:21 PST 2016
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>
> email: 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>> 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>
>
> 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.887
> -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>
> email: 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>
> 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