Anelastic liquid approximation

[antonio-schettino]

Hi,

I am trying to implement a system based on the anelastic liquid approximation by the following calls:

divergence = -alpha * g * velocityField[1] * S_L / (Gamma * cp)
advDiff = uw.systems.AdvectionDiffusion(phiField       = temperatureField,
                                        velocityField  = velocityField,
                                        fn_diffusivity = diffusivity,
                                        conditions     = tempBC,
                                        method         = 'SLCN')
stokes = uw.systems.Stokes(velocityField = velocityField,
                           pressureField = pressureField,
                           conditions    = velocityBC,
                           fn_viscosity  = viscosityFn,
                           fn_bodyforce  = buoyancyFn,
                           fn_source     = divergence)
solver = uw.systems.Solver(stokes)
solver.set_inner_method("lu")
solver.solve(nonLinearIterate = True,nonLinearMaxIterations = 20,nonLinearTolerance = 0.02)

where S_L is the distance scaling factor. In this instance, I have no bulk viscosity. However, the resulting velocity field has a divergence field in a range two/three orders of magnitude less than the expected values stored in the variable 'divergence':

gradV = velocityField.fn_gradient
div = gradV[0] + gradV[3]
maxdiv1 = np.max(div.evaluate(mesh)) * S_v_cmyr * 1e3 / S_L
mindiv1 = np.min(div.evaluate(mesh)) * S_v_cmyr * 1e3 / S_L
maxdiv2 = np.max(divergence.evaluate(mesh)) * S_v_cmyr * 1e3 / S_L
mindiv2 = np.min(divergence.evaluate(mesh)) * S_v_cmyr * 1e3 / S_L
print(mindiv1,maxdiv1)
print(mindiv2,maxdiv2)

Result:

-0.0120159773895 0.0108061217017
-0.000133182970101 3.66106035129e-05

Is there anybody who can help me to find the mistake?

[jmansour]

@antonio-schettino do you get the expected result for a constant divergence?

[lmoresi]

Julian Giordani has benchmarked these cases and should be able to discuss this with you.

Prof Louis Moresi

louis.moresi@anu.edu.au

(m) +61 4 0333 1413

(us) +1 505 349 4425

www.moresi.infohttp://www.moresi.info/

www.facebook.com/underworldcodehttp://www.facebook.com/underworldcode

@LouisMoresihttps://twitter.com/LouisMoresi

On 12 Nov 2019, 17:12 +1100, John Mansour notifications@github.com, wrote:

@antonio-schettinohttps://github.com/antonio-schettino do you get the expected result for a constant divergence?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/underworldcode/underworld2/issues/437?email_source=notifications&email_token=ADABPI3DHBZHFKPC66OA4XTQTJCLZA5CNFSM4JLVCEK2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEDZE3PI#issuecomment-552750525, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADABPIZO5LUA6622CPEJMJ3QTJCLZANCNFSM4JLVCEKQ.

[antonio-schettino]

Assuming constant vz(x,y) = 1 cm/yr in the source term, we would have divergence = -54.36677777777778 (in nondimensional units). Using such constant source term, I obtain an effective divergence ranging between -3566.80400111 and 2627.01917769 (in nondimensional units). Here I am assuming Gamma = 1 (Gruneisen’s parameter), alpha = 3.5e-5 K-1, cp = 1200.0 J K-1 kg-1, S_L = 670.0e3 m.

When vz(x,y) = 0 cm/yr, the resulting divergence field varies between -2287.15588148 and 2034.92352961. The distribution seems related to the vertical component of the velocity field anyway.

[antonio-schettino]

I have not solved even introducing the parameter lambda. The following call:

stokes = uw.systems.Stokes(velocityField = velocityField, pressureField = pressureField, conditions = velocityBC, fn_viscosity = viscosityFn, fn_bodyforce = buoyancyFn, fn_source = divergence, fn_one_on_lambda = fn.misc.constant(1.0 / 1e22))

produces anyway an incorrect result:

maxdiv1 = np.max(div.evaluate(mesh)) mindiv1 = np.min(div.evaluate(mesh)) maxdiv2 = np.max(divergence.evaluate(mesh)) mindiv2 = np.min(divergence.evaluate(mesh)) print(mindiv1,maxdiv1) print(mindiv2,maxdiv2)

-2283.4976633 2052.74961344 -25.309347956 6.96460440024

[julesghub]

The Anelastic approximation (ALA) was not tested. We only tested the TALA approximation for compressible convection models. These tests were made with an older version of the code and aren't currently in our testing suit.

There is a simple test for the use of fn_one_on_lambda at https://github.com/underworldcode/underworld2/blob/master/docs/test/compress.ipynb

@antonio-schettino, interesting that vz=0 produces variations. That is likely because divergence is a function of vz by definition - this will activate the non linear picard iterations . Can you set fn_source to a constant, removing the dependence on vz. Does that help?

[jmansour]

Also, what is your boundary condition? Note that if you've added a source, you must also allow an open/outflow bc or some other compliance in the system (such as a compressible region) otherwise you will violate volume conservation.

Edit: Note that this is if you have a net positive/negative source. Note also that while you can prescribe a balance of in/outflows and sources/sinks, this can be fragile due to the discrete nature by which these configurations are captured. In practise it is usually more robust to allow compliance in the system via an amount of compressibility somewhere or an open boundary.

[antonio-schettino]

In the original model, there is outflow and inflow, but I check at each step that the net flux is zero. Furthermore, I tested the model with free slip BCs as well. For example, this is the result with imposed vx = const at the top boundary:

In the case of inflow/outflow along the right boundary, we have: Finally, assuming a constant source term (vy = 0.1 cm/yr everywhere):

presudovisc      = 1.0e24
fn_one_on_lambda = fn.misc.constant(S_eta / presudovisc)
source           = -alpha * g * 0.1 * S_L / (Gamma * cp * S_v_cmyr) + pressureField * fn_one_on_lambda 
stokes = uw.systems.Stokes(velocityField     = velocityField,
                           pressureField     = pressureField,
                           conditions        = [velocityBC,nbc_bottom],
                           fn_viscosity      = viscosityFn,
                           fn_bodyforce      = buoyancyFn,
                           fn_source         = source,
                           fn_one_on_lambda  = fn_one_on_lambda)

the result does not change significantly:

gradV = velocityField.fn_gradient
div1 = gradV[0] + gradV[3]
div2 = -alpha * g * fn.misc.constant(0.1) * S_L / (Gamma * cp * S_v_cmyr)
maxdiv1 = np.max(div1.evaluate(mesh))
mindiv1 = np.min(div1.evaluate(mesh))
maxdiv2 = np.max(div2.evaluate(mesh))
mindiv2 = np.min(div2.evaluate(mesh))
print(mindiv1,maxdiv1)
print(mindiv2,maxdiv2)

-2275.85615183 1991.14071067
-5.43667777778 -5.43667777778

[antonio-schettino]

Another interesting issue: Setting H = 0 (the source term in Eq. 1.1 of the Underworld documentation) in the compressible equation (see the test suggested by Julian at https://github.com/underworldcode/underworld2/blob/master/docs/test/compress.ipynb) the divergence should satisfy: div v + p/Lambda = 0 everywhere. Instead, we have:

So I guess, the problem is the pressure? I have initialized the density using the Adams-Williamson equation and T adiabatic plus axial anomaly. The resulting pressure field is:

Maybe the correct pressure in the equation div v + p/Lambda = 0 should be the total pressure?

[jmansour]

With your last results, I'll note that the results look like the divergence is pretty much zero everywhere, except for the noise around the centre. Keep in mind that the numerical divergence is quite difficult to calculate accurately, as it relies on the velocity gradients which are discontinuous and therefore quite noisy. As such the artifacts in the centre may or may not be of concern (as far as the divergence is concerned).

A better approximation may possibly be found using the mesh centroids. What are the results for the following:

maxdiv1 = np.max(div1.evaluate(mesh.subMesh))
mindiv1 = np.min(div1.evaluate(mesh.subMesh))
maxdiv2 = np.max(div2.evaluate(mesh.subMesh))
mindiv2 = np.min(div2.evaluate(mesh.subMesh))
print(mindiv1,maxdiv1)
print(mindiv2,maxdiv2)

here I'm using evaluate(mesh.subMesh) which will instead sample the fields at element centroids, as opposed to evaluate(mesh) which will sample fields at element edges.

Can you also report your solver iteration counts?

[antonio-schettino]

John, this is a good point. Below are listed the results in the case of my anelastic liquid approximation (ALA) test, where the theoretical divergence is calculated on the basis of the vertical component of the velocity:

presudovisc      = fn.misc.constant(S_eta / 1.0e25)
source           = -2.0 * alpha * g * velocityField[1] * presudovisc * viscosityFn * S_L / (3.0 * Gamma * cp )
stokes = uw.systems.Stokes(velocityField     = velocityField,
                           pressureField     = pressureField,
                           conditions        = velocityBC,
                           fn_viscosity      = viscosityFn,
                           fn_bodyforce      = buoyancyFn,
                           fn_one_on_lambda  = presudovisc,
                           fn_source         = source)

As you can see, the difference between solver and theoretical divergences is two orders of magnitude.

The counts are:

div1 = gradV[0] + gradV[3]
div2 = alpha * g * velocityField[1] * S_L / (Gamma * cp)
maxdiv1 = np.max(div1.evaluate(mesh.subMesh))
mindiv1 = np.min(div1.evaluate(mesh.subMesh))
maxdiv2 = np.max(div2.evaluate(mesh.subMesh))
mindiv2 = np.min(div2.evaluate(mesh.subMesh))
print(mindiv1,maxdiv1)
print(mindiv2,maxdiv2)

-630.25989929 526.852286167
-6.21654775112 15.1558698752

If we eliminate the source term, so that div v + p/lambda = 0, the results are similar:

Counts:

-717.686954814 657.981846989
-9.9914125649 108.335458823

My point is: the pressure used in these calculations is the dynamic pressure stored in the system variable pressureField. However, in the formulation of Hughes (1987) the variable p in the constitutive equation (and in the equation div v + p/lambda = 0) is the total pressure. Is this the error source? If I estimate div v by (p + p_ref)/lambda, the discrepancy is reduced by one order of magnitude:

Counts:

-630.409559082 527.038562033
-1404.86245805 95.8208239534

[antonio-schettino]

I have fixed a bug in the scaling factors, now the solver works fine. Unfortunately, it seems there is no way to implement ALA in underworld.

[jmansour]

I don't see any fundamental reason why ALA cannot be implemented, but it's hard to know where your model is going wrong.

If you wish to pursue it, I'd suggest you simplify everything as much as possible, including numbers close to 1 to ensure this isn't a problem of sensitivity.

Here, I created a simple test based on our RT example. I've called it rt_ala.ipynb, but really it simply just sets fn_source=velocityField[1].

https://mybinder.org/v2/gist/jmansour/e0ddd91590b7152efd13034405994c15/master

It's open at the top to accommodate the divergence. I've used "Q2/dpc1" elements as these seemed to give better results, although "Q1/dq0" wasn't terrible.

Note that the sign of the source is actually reversed with respect to the documentation.

[antonio-schettino]

The problem, in my opinion and assuming that I have understood the documentation, is that underworld does not consider the compressional term (divergence of v) in the costitutive equation:

where I have assumed that Stokes condition holds. I other words, I understand that the stress is calculated assuming div v = 0. Consequently, the penalty trick (div v + p/lambda = H) can only deal with quasi-incompressible conditions. Is this interpretation correct?

[jmansour]

Yes @antonio-schettino this is exactly correct.

[lmoresi]

Just a short comment on this. The following paper does make use of non-zero divergence term and there are quite a few issues when it comes to non-linearity when the divergence is implemented as an explicit source term.

Muhlhaus, H.-B., and L. N. Moresi. “A Model for Trap Door Flow from a Deep Container.” In Desiderata Geotechnica, edited by Wei Wu, 113–18. Springer Series in Geomechanics and Geoengineering. Springer International Publishing, 2019.

In the original formulation of this algorithm (Moresi et al, 2002/3) , we did implement compressible visco-elasticity and hence a second viscosity term of the form required. The constitutive behaviour in that case uses the same matrices as the penalty term does. This has not been implemented in the current underworld.

Moresi, L., F. Dufour, and H.-B. Mühlhaus. “A Lagrangian Integration Point Finite Element Method for Large Deformation Modeling of Viscoelastic Geomaterials.” Journal of Computational Physics 184, no. 2 (2003): 476–97. doi.org/10.1016/S0021-9991(02)00031-1.

———. “Mantle Convection Modeling with Viscoelastic/Brittle Lithosphere:” Pure and Applied Geophysics 159, no. 10 (2002): 2335–56. doi.org/10.1007/s00024-002-8738-3.

On 2 Dec 2019, 8:30 AM +0000, antonio-schettino notifications@github.com, wrote:

Closed #437https://github.com/underworldcode/underworld2/issues/437.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/underworldcode/underworld2/issues/437?email_source=notifications&email_token=ADABPI6OY26LU6M3ECTHJW3QWTBP7A5CNFSM4JLVCEK2YY3PNVWWK3TUL52HS4DFWZEXG43VMVCXMZLOORHG65DJMZUWGYLUNFXW5KTDN5WW2ZLOORPWSZGOVGHNRTQ#event-2844711118, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADABPI3AAA3IUP3WD7PTWFLQWTBP7ANCNFSM4JLVCEKQ.