[CIG-SHORT] Stress issue with spherical model

[Demian Gomez]

Brad,

How do I "verify what initial stresses PyLith is using by turning off
gravity and not having any deformation"? I can turn off gravity and let the
displacement at the fault be = 0, but how do I constrain the displacements
everywhere in the model so that after I run Pylith the stresses are = to
the initial stresses?

Also, the idea of querying the spatialdb with a Python script is good, but
unfortunately, the example you refer to (exodus_add_properties.py) it's too
simple for my case. After several hours of trial and error I wasn't able to
query the spatialgriddb using the coordinates of my mesh. I think the
problem is related to the coordinate system. From the example, I've defined
my coordinate system object as:

    # Coordinate system for mesh (must match coordsys in ExodusII file)
    cs = CSGeo()
    cs.inventory.ellipsoid = "sphere"
    cs.inventory.datum_vert = "ellipsoid"
    cs.inventory.is_geocentric = True
    cs._configure()
    cs.initialize()

When I try to query the geocentric coordinates of the model, I do not get
any values from the db (vector err full of ones). If I use points in lat
lon h, then the query returns values (err vector full of zeros). I've
attached the python script based on the example you suggested.

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 28, 2016 at 1:17 PM, Brad Aagaard <baagaard at usgs.gov> wrote:

> 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
>>
>>
> _______________________________________________
> 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/20161128/1dfbff3c/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: query_db.py
Type: text/x-python
Size: 2497 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20161128/1dfbff3c/attachment.py>