[CIG-SHORT] nondimensionalization question

Tue Aug 21 10:57:06 PDT 2012

Hello Brad,

Thanks a lot for your help. I just implemented the fix you proposed and I can confirm that everything behaves as expected now. 

Regards,
Lucas
________________________________________
From: Brad Aagaard [baagaard at usgs.gov]
Sent: Tuesday, August 21, 2012 13:11
To: Lucas Abraham Willemsen; cig-short at geodynamics.org
Subject: Re: [CIG-SHORT] nondimensionalization question

Lucas,

This issue turned out to be a subtle bug related to using
NondimElasticQuasistatic with explicit time stepping. Thanks for bring
this issue to our attention and pushing us to look into it. This helps
make PyLith more robust for everyone. I have updated
NondimElasticQuasistatic as described below. This fix will be included
in the next release. You can manually patch your installed version
and/or binary as described below.

If we start with the NondimElasticDynamic scales, we have:

length_scale = vs * Tmin
pressure_scale = density * vs**2 = shear_modulus
time_scale = Tmin
density_scale = density

For NondimElasticQuasistatic, we have:

length_scale = L
pressure_scale = shear_modulus
time_scale = relaxation_time

For the scales to match, we have:

L = vs * Tmin
shear_modulus = density * vs**2
relaxation_time = Tmin

So clearly, the density_scale in the NondimElasticQuasistatic should be
density = shear_modulus / vs**2, where vs = L / relaxation_time. Using a
default value of 1.0 will result in incorrect nondimensionalization of
density values.

I have added setting the density scale to shear_modulus / vs**2 in
NondimElasticQuasistatic. You can manually change
NondimElasticQuasistatic.py in the installed location to match:

http://www.geodynamics.org/wsvn/cig/cs/spatialdata/trunk/spatialdata/units/NondimElasticQuasistatic.py?op=diff&rev=0&sc=0

I added the following lines

     # Compute implicit density scale.
     vs = self.inventory.lengthScale / self.inventory.relaxationTime
     self.setDensityScale(self.inventory.shearModulus / vs**2)

Regards,
Brad

On 08/21/2012 09:29 AM, Lucas Abraham Willemsen wrote:
> Hi Brad,
>
> In the project I linked (http://web.mit.edu/lawillem/www/nondimtest.zip) I use an explicit timescale. It uses a spontaneous rupture with slip weakening, the initial_traction on the fault is slightly larger than the strength (76.001 MPa vs 76.000 MPa). The same issue can be observed by using other frictional models. No matter what values I choose for the dynamic object I always get a displacement_x of 0.0424 meter at time 0.495 seconds. In this specific case I will use the following dynamic values:
>
> shear_wave_speed = 3000*m/s
> mass_density = 1.0*kg/m**3
> wave_period = 0.50*s
>
> I chose a mass density of 1.0*kg/m**3 since the quasistatic object defaults to this. Any other mass_density value (say 2500) would still give the same displacement_x. These dynamic nondim values corresponds to the following scales:
>
> length = vs * period = 3000 * 0.5 = 1500
> time = 0.5
> pressure = vs**2 * density = 3000**2 * 1.0 = 9000000
> density = 1.0
>
> I then uncomment the quasistatic object in the linked project (and comment out the dynamic). I use the corresponding values:
>
> ----------------------------------
> length_scale = 1500*m
> shear_modulus = 9000000*Pa
> relaxation_time = 0.5*s
> ----------------------------------
>
> And the final displacement_x is still 0.0424m as expected. Like you also mentioned, it is possible to get the same results using both objects. But unlike the dynamic object, any change in the ratio between the values in the quasistatic object will alter the final solution dramatically. I will demonstrate by changing only the pressure scale, and by only 1%.
>
> ----------------------------------
> MODIFIED VALUES
> ----------------------------------
> length_scale = 1500*m
> shear_modulus = 9100000*Pa
> relaxation_time = 0.5*s
> ----------------------------------
>
> The displacement_x after 0.495 seconds is now 0.0404m. A change of almost 5% compared to the previous 0.0424m. My point is that it is possible to get the 'correct' (not varying) behavior of the dynamic object using a quasistatic object, but even the slightest deviation from these values will give wrong results. My intuition says that a change of only 1% in for instance the pressure scale should not change the final result since the dimensionless values have the same order of magnitude. After the simulation is complete they will be made dimensional again and the result should be the same.
>
> (if I would have changed the length_scale or the relaxation_time instead of the shear_modulus, I would have observed the same behavior)
>
> I hope this example helps.
> Thanks for your help,
> Lucas
>
> ________________________________________
> From: Brad Aagaard [baagaard at usgs.gov]
> Sent: Tuesday, August 21, 2012 11:35
> To: Lucas Abraham Willemsen
> Cc: cig-short at geodynamics.org
> Subject: Re: [CIG-SHORT] nondimensionalization question
>
> Lucas,
>
> In general, we setup NondimElasticQuasistatic so that the scales are
> more user-friendly for quasi-static problems using implicit time
> stepping. Likewise, the scales in NondimElasticDynamic are more
> user-friendly for dynamic problems with explicit time stepping. You can
> use whichever is more convenient. Choosing the proper scales for a given
> problem is an independent issue.
>
> Here are some issues to keep in mind related to nondimensionalization
> and PyLith:
>
> 1. The implicit formulation is much more sensitive to the nondimensional
> scales because it affects the conditioning of the linear system and the
> absolute tolerances in the iterative solve. The solution vector contains
> both displacements and the Lagrange multipliers (fault tractions), which
> ideally will be the same order of magnitude when nondimensionalized. If
> they differ by more than a few orders of magnitude, the system will be
> very ill-conditioned and convergence will be slow and with improper
> solver tolerances the errors will be large.
>
> 2. The nondimensional scales affect the tolerances associated with
> determining whether a small amount of slip is above or below the zero
> threshold. Improper selection of the tolerances and scales can lead to
> slip when the fault should stick and vice versa.
>
> 3. If the *only* change in a simulation is between
> NondimElasticQuasistatic and NondimElasticDynamic and *all* of the
> scales match, the solutions should be identical. If the solutions
> differ, then there is a bug. I will check this today.
>
>
> Can you be more precise in what parameters you are changing? In the test
> cases you are running, are you always using explicit time stepping? Are
> you using spontaneous rupture or prescribed slip?
>
> Brad
>
> On 08/21/2012 07:48 AM, Lucas Abraham Willemsen wrote:
>> Hi Brad,
>>
>> I just tested whether the numerical damping affected the results, but
>> it turns out it does not in this case. The dynamic nondim object
>> consistently gives the same results independent of the values chosen,
>> while the results with the quasistatic nondim object differ depending
>> on the values chosen (up to more than 10 orders of magnitude).
>>
>> Instead of spending more time on this it may be more fruitful if I
>> just use the dynamic object for both my dynamic (explicit) and
>> quasistatic (implicit) formulations, since it seems to consistently
>> give the same result independent of chosen values, .
>>
>> thanks, Lucas ________________________________________ From: Brad
>> Aagaard [baagaard at usgs.gov] Sent: Monday, August 20, 2012 23:55 To:
>> Lucas Abraham Willemsen Cc: cig-short at geodynamics.org Subject: Re:
>> [CIG-SHORT] nondimensionalization question
>>
>> Lucas,
>>
>> You can turn off the numerical damping by setting it to a very small
>> positive number (it looks like we incorrectly disallow exactly 0.0):
>>
>> formulation = pylith.problems.ExplicitLumped
>> formulation.norm_viscosity = 1.0e-30
>>
>> Brad
>>
>>
>> On 8/20/12 8:46 PM, Lucas Abraham Willemsen wrote:
>>> Hi Brad,
>>>
>>> I realize that in the example I provided the dynamic and the
>>> quasi-static scaling parameters are different. It is possible to
>>> get the same results using both objects.
>>>
>>> It is indeed possible to mimic the results of the dynamic object
>>> using a quasi-static object. But my problems is as follows:
>>> Changing any of the scaling values individually in the dynamic
>>> object will yield the same result. Changing any of the scaling
>>> values individually in the quasi-static object will change the
>>> final result. I don't understand why that happens. Could it be the
>>> nondimensionalized numerical damping you mentioned?
>>>
>>> best, Lucas ________________________________________ From:
>>> cig-short-bounces at geodynamics.org
>>> [cig-short-bounces at geodynamics.org] on behalf of Brad Aagaard
>>> [baagaard at usgs.gov] Sent: Monday, August 20, 2012 23:35 To:
>>> cig-short at geodynamics.org Subject: Re: [CIG-SHORT]
>>> nondimensionalization question
>>>
>>> Lucas,
>>>
>>> Here are the expressions for the length scales for the
>>> nondimensionalization:
>>>
>>> NondimElasticDynamic
>>>
>>> length_scale = vs * period pressure_scale = vs**2 * density
>>> time_scale = period density_scale = density
>>>
>>> NondimElasticQuasistatic
>>>
>>> length_scale = lengthScale pressure_scale = shearModulus time_scale
>>> = relaxationTime
>>>
>>> In your pylithapp.cfg I see
>>>
>>> normalizer = spatialdata.units.NondimElasticDynamic
>>> [pylithapp.timedependent.normalizer] shear_wave_speed = 5200*m/s
>>> mass_density = 1.9e+3*kg/m**3 wave_period = 0.34*s
>>>
>>> #normalizer = spatialdata.units.NondimElasticQuasistatic
>>> #[pylithapp.timedependent.normalizer] #length_scale = 1.2*m
>>> #shear_modulus = 1.0e+10*Pa #relaxation_time = 1.3*s
>>>
>>> Note that the NondimElasticQuasistatic object does not provide a
>>> density scale (because inertia does not enter into quasi-static
>>> problems). The default density_scale is 1.0.
>>>
>>> Using the above expressions, your length_scale, pressure_scale,
>>> time_scale, and density_scale do not match. The length_scale,
>>> pressure_scale, and time_scale differ by less than an order of
>>> magnitude, but the density scales obviously differ by 1.9e+3. This
>>> almost surely explains the difference in results.
>>>
>>> Additionally, in the explicit time stepping we include
>>> nondimensional numerical damping, which will vary if the time
>>> scales are different.
>>>
>>> Regards, Brad
>>>
>>> On 8/20/12 7:17 PM, Lucas Abraham Willemsen wrote:
>>>> Hi Charles,
>>>>
>>>> The layout of the files you linked is different, but it seems
>>>> like they apply the scales in the same way. The behavior of the
>>>> quasistatic nondimensionalization object seems to be strange, at
>>>> least when a dynamic problem is used. The test case I presented
>>>> was using Pylith 1.7.1 by the way.
>>>>
>>>> best, Lucas
>>>> ------------------------------------------------------------------------
>>>>
>>>>
> *From:* Charles Williams [willic3 at gmail.com]
>>>> *Sent:* Monday, August 20, 2012 22:07 *To:* Matthew Knepley *Cc:*
>>>> Lucas Abraham Willemsen; cig-short at geodynamics.org *Subject:* Re:
>>>> [CIG-SHORT] nondimensionalization question
>>>>
>>>> Also, try looking in the trunk, which is the current
>>>> implementation. I'm not sure how much might have changed:
>>>>
>>>> http://geodynamics.org/svn/cig/cs/spatialdata/trunk/spatialdata/units/
>>>>
>>>>
>>>>
> Cheers,
>>>> Charles
>>>>
>>>>
>>>> On 21/08/2012, at 1:21 PM, Matthew Knepley wrote:
>>>>
>>>>> On Mon, Aug 20, 2012 at 7:19 PM, Lucas Abraham Willemsen
>>>>> <lawillem at mit.edu <mailto:lawillem at mit.edu>> wrote:
>>>>>
>>>>> Hi,
>>>>>
>>>>> In preparation for my research I am currently investigating
>>>>> the nondimensionalization implementation in Pylith. The way I
>>>>> understand things is that nondimensionalization should have no
>>>>> effect as long as the double precision is good enough in
>>>>> preventing round-off errors. I used a dynamic simulation with
>>>>> a
>>>>>
>>>>>
>>>>> This is not the whole story. First, you are not just competing
>>>>> against roundoff error, but also truncation error, and these
>>>>> errors can be inflated by the condition number of your problem,
>>>>> so poor scaling can result in very wrong answers. Second, poor
>>>>> scaling can result in very poor solver performance as well.
>>>>>
>>>>> slip-weakening fault in order to test this theory. When I use
>>>>> a dynamic nondimensionalization object everything is as
>>>>> expected. Changing the scales 'shear_wave_speed',
>>>>> 'mass_density' and 'wave_period' has no effect on the final
>>>>> displacement.
>>>>>
>>>>> I was told (maybe in error?) that both dynamic and quasistatic
>>>>> nondimensionazation objects are valid in a dynamic simulation
>>>>> and should give the same result (they do the same thing?). But
>>>>> when I test this it does not work. changing any of the scales
>>>>> in the quasistatic nondimensionalization object changes the
>>>>> final displacements significantly (orders of magnitude).
>>>>>
>>>>>
>>>>> To make sure that you are doing what you want, you must switch
>>>>> to a direct solver. In this case, it means using FieldSplit,
>>>>> full Schur factorization, LU for the displacements, and a very
>>>>> low tolerance (1e-12 or so) for the fault tractions. We should
>>>>> definitely make an options file for these choices. With these,
>>>>> we will know whether solver convergence is influencing your
>>>>> results.
>>>>>
>>>>> Matt
>>>>>
>>>>> I browsed around in the source code for the implementations of
>>>>> the objects and I found this. It does seem like both these
>>>>> objects essentially do the same thing (except for the fact that
>>>>> the quasistatic object will always have a default density scale
>>>>> since it is intended for quasistatic problems).
>>>>>
>>>>> http://geodynamics.org/svn/cig/cs/spatialdata/tags/v0.5.2/spatialdata/units/NondimElasticDynamic.py
>>>>>
>>>>>
> http://geodynamics.org/svn/cig/cs/spatialdata/tags/v0.5.2/spatialdata/units/NondimElasticQuasistatic.py
>>>>>
>>>>> (it says v0.5.2 in the link, but there is no higher one. Is
>>>>> this what is currently used? Could not find the files in the
>>>>> source code download)
>>>>>
>>>>> A test project can be found here:
>>>>> http://web.mit.edu/lawillem/www/nondimtest.zip
>>>>>
>>>>> Note how changing the dynamic scales changes nothing, while
>>>>> the quasistatic ones do influence the final displacements
>>>>> significantly.
>>>>>
>>>>> best, Lucas
>>>>>
>>>>> P.S. My motivation for this question is that I plan to
>>>>> investigate the difference between a quasi-static and dynamic
>>>>> simulation with rate and state friction. In order to make the
>>>>> transition from dynamic to real quasistatic (implicit
>>>>> formulation) I first wanted to change the nondimensionalization
>>>>> object to quasistatic (while problem remains dynamic, explicit
>>>>> timestep) and get the same results.
>>>>>
>>>>> _______________________________________________ CIG-SHORT
>>>>> mailing list CIG-SHORT at geodynamics.org
>>>>> <mailto:CIG-SHORT at geodynamics.org>
>>>>> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> -- What most experimenters take for granted before they begin
>>>>> their experiments is infinitely more interesting than any
>>>>> results to which their experiments lead. -- Norbert Wiener
>>>>> _______________________________________________ CIG-SHORT
>>>>> mailing list CIG-SHORT at geodynamics.org
>>>>> <mailto:CIG-SHORT at geodynamics.org>
>>>>> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>>>>
>>>> Charles A. Williams Scientist GNS Science 1 Fairway Drive,
>>>> Avalon PO Box 30368 Lower Hutt 5040 New Zealand ph (office):
>>>> 0064-4570-4566 fax (office): 0064-4570-4600 C.Williams at gns.cri.nz
>>>> <mailto:C.Williams at gns.cri.nz>
>>>>
>>>>
>>>>
>>>> _______________________________________________ CIG-SHORT mailing
>>>> list CIG-SHORT at geodynamics.org
>>>> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>>>>
>>>
>>> _______________________________________________ CIG-SHORT mailing
>>> list CIG-SHORT at geodynamics.org
>>> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>>>
>>
>>
>
>