[CIG-SHORT] Stress issue with spherical model

[Brad Aagaard]

Demian,

I don't quite follow your through process in trying to verify your 
problem setup. The process I would use is (1) verify that the gravity 
field you are getting is what you expect, (2) the initial stresses you 
are specifying in the spatial database match what PyLith is using, (3) 
check to make sure the initial stresses balance the gravitational body 
forces.

For (1) un a geocentric coordinate system, GravityField will use the 
radial direction from the origin as the direction of the gravity field. 
This needs to match your assumptions for setting up the initial stresses.

For (2) as Charles mentioned you can verify what initial stresses PyLith 
is using by turning off gravity and not having any deformation. The 
stress field output should match your initial stresses as given in the 
spatial database.

If you have questions about the transformations among coordinate 
systems, use the Proj.4 command line tools cs2cs and proj.

You can also use the Python interface to run queries of the 
GravityField, do coordinate transformations, etc. See 
examples/meshing/exodus_add_properties.py for an example of doing a 
spatial database query from a Python script.

Regards,
Brad


On 11/24/2016 08:40 AM, Demian Gomez wrote:
> Dear all,
>
> I think I found the problem with the lithostatic stress in my spherical
> model. I'm not sure if it's a bug or not. I am specifying the stress
> using a SpatialGridDB in lat lon h, where I declare my coordinate system as:
>
>   cs-data = geographic {
>   to-meters = 1
>   space-dim = 3
>   ellipsoid = sphere
>   datum-vert = ellipsoid
>   is-geocentric = false
> }
>
> In pylithapp.cfg I declare my coordinate system as:
>
> [pylithapp.mesh_generator.reader]
> coordsys = spatialdata.geocoords.CSGeo
> coordsys.ellipsoid = sphere
> coordsys.datum_vert = ellipsoid
> coordsys.is_geocentric = true
>
> The problem seems to be that Pylith is still using an ellipsoid instead
> of using ellipsoid = sphere for the SpatialGridDB (I think). I verified
> this by doing the following two experiments (see attached PDF):
>
> 1) I've created the spatialdb by calculating the lithostatic stress on a
> lat lon h grid.
> 2) Same as before but this time I added a stress "offset" of 10 m
> (offset = 10*g*rho). This is equivalent to radially scaling the spatialdb.
>
> From the attached figure, you'll notice that in experiment 1, the
> displacements are 0 at the center of the domain (the picture shows the
> domain cut by a plane). In experiment 2, which has the stress offset,
> the center of the domain shows displacements pointing towards the
> geocenter (stress defect) and a zero displacement ring around the center
> of the domain. Moving away from this ring (towards the edge of the
> domain), there is a region where the displacements point up, showing an
> excess of stress. I've depicted each experiment schematically in the figure.
>
> Would it be possible to verify if what I found is correct? Or maybe
> there's something wrong with my configuration.
>
> Thanks!
> 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
>