[CIG-SHORT] Spherical 3d gravity initial stress

Tue Nov 8 16:08:25 PST 2016

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
email: gomez.124 at osu.edu

On Mon, Nov 7, 2016 at 12:42 PM, Brad Aagaard <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
>>
>> 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.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::vect
>> or<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
>> email: gomez.124 at osu.edu <mailto:gomez.124 at osu.edu>
>>
>>
>>
>> _______________________________________________
>> CIG-SHORT mailing list
>> CIG-SHORT at geodynamics.org
>> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20161108/c792c848/attachment-0001.html>