[CIG-SHORT] Spherical 3d gravity initial stress

[Demian Gomez]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20161106/5685e644/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: geometry_displacement.png
Type: image/png
Size: 1611482 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20161106/5685e644/attachment-0001.png>