[CIG-SHORT] Spherical 3d gravity initial stress

[Brad Aagaard]

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::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
> 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
>