[CIG-CS] Discontinuous viscosity with curved interfaces
Walter Landry
walter at geodynamics.org
Tue May 8 11:02:43 PDT 2012
Hello Everyone,
Back in March, I posted some results for a Stokes solver for
discontinuous viscosities.
http://www.geodynamics.org/pipermail/cig-cs/2012-March/000060.html
It turns out that the method I implemented was second order convergent
in pressure. That is more than is needed, especially for adaptive
codes. So I have simplified the method to first order in pressure,
second order in velocity. This also shrinks the stencil which should
help a parallel implementation. I think the current implementation
only really needs 2 ghost points, instead of the traditional 1 ghost
point.
The rate of convergence seems unaffected. Once again, I ran the SOLCX
benchmark [1] on a 64x64 grid with the jump at x=0.4
Viscosity ratio # iterations
1 18132
1e2 22069
1e4 22142
1e6 22142
1e8 22142
1e10 22142
The number of iterations scales as the number of grid points, as
expected for Gauss-Seidel. For a jump of 10^10
h # iterations
1/32 5443
1/64 22142
1/128 89595
The L_Infinity, L1, and L2 errors in the pressure and velocity for a
jump of 10^10 are
p
h L_inf(p) L2(p) L1(p)
1/32 0.15296 0.0481187 0.0336702
1/64 0.08159 0.0249955 0.0169148
1/128 0.03008 0.00913483 0.0064351
vx
h L_inf(vx) L2(vx) L1(vx)
1/32 0.000217779 3.59084e-05 1.47698e-05
1/64 7.32885e-05 1.1595e-05 5.16648e-06
1/128 0.00014967 4.93944e-05 2.54406e-05
vy
h L_inf(vy) L2(vy) L1(vy)
1/32 0.00132709 0.00025398 0.000103956
1/64 0.000418997 7.61899e-05 3.17731e-05
1/128 0.0002025 6.3144e-05 3.1889e-05
The non-convergent L_inf(vx) is because the location of the interface
relative to the grid is moving around. It can be closer to the high
viscosity grid point or the low viscosity grid point. The method
converges regardless, but with different pre-factors. I can verify
this by setting the jump to be at x=(13-1/16)/32.
p
h L_inf(p) L2(p) L1(p)
1/32 0.15372 0.0483512 0.0337975
1/64 0.08315 0.0254631 0.0171066
1/128 0.04295 0.0129933 0.00856586
vx
h L_inf(vx) L2(vx) L1(vx)
1/32 0.000219854 3.6038e-05 1.46592e-05
1/64 8.97136e-05 1.37678e-05 4.8201e-06
1/128 2.61837e-05 3.90678e-06 1.37179e-06
vy
h L_inf(vy) L2(vy) L1(vy)
1/32 0.00131613 0.000251982 0.00010304
1/64 0.00037517 6.97169e-05 2.81431e-05
1/128 0.000102199 1.85925e-05 7.52455e-06
I have also generalized the code for arbitrary curved interfaces. To
test it, I solved the problem of a 2D circular inclusion in pure
shear [2]. The inclusion has a high viscosity and the medium is low
viscosity. For a 64x64 grid, with an inclusion r_c=0.2 and the
outer boundary at 1, the number of iterations does not change
significantly.
Viscosity # Iterations
Ratio
1 25287
1e2 23310
1e4 23195
1e6 23195
1e8 23195
1e10 23195
It is interesting that the number of iterations actually decreases as
viscosity increases. With small viscosity, the solution is just a
linear ramp. Gauss-Seidel is not effective at smoothing this pure
long-wavelength solution. Multigrid should remedy this.
For a viscosity ratio of 1e10, the L_infinity, L2, and L1 errors are
p
h L_inf(p) L2(p) L1(p)
1/32 1.40381 0.0846449 0.0469965
1/64 1.046 0.0748157 0.0672745
1/128 0.659582 0.0453235 0.040746
vx
h L_inf(vx) L2(vx) L1(vx)
1/32 0.0074005 0.000856215 0.000510485
1/64 0.00362048 0.000609875 0.000437284
1/128 0.002967 0.000504043 0.000346278
vy
h L_inf(vy) L2(vy) L1(vy)
1/32 0.0153628 0.00111183 0.000576541
1/64 0.0047258 0.000615105 0.000439425
1/128 0.002967 0.000505353 0.000347465
It is a bit difficult to interpret these numbers because at the higher
resolutions, the error in the velocity becomes dominated by errors in
the outer boundary. Phase errors may also be playing a role, since I
am using level sets to determine where the interface is. I am
attaching a picture of the error in pressure for all three
resolutions. I am also attaching side-by-side pictures of the
pressure and the analytic solution for all three resolutions (left is
numerical, right is analytic).
So this is working pretty well. I do not see any real problems with
it. The next step is to put this into Gamr, my parallel, adaptive
code. Then I can apply multigrid and see how that fares. I have to
fix the coarsening and refinement operators to handle the jumps, but
that appears to be straightforward. Then it would be possible to do
an apples-to-apples comparison with Aspect.
Cheers,
Walter Landry
[1] Analytic solutions for Stokes$B!G(B flow with lateral variations in viscosity
Shijie Zhong
Geophys. J. Int.(1996) 124, 18-28
Section 2.3
[2] Analytical solutions for deformable elliptical inclusions
in general shear
Daniel W. Schmid and Yuri Yu. Podladchikov
Geophys. J. Int. (2003) 155, 269$(Q#|(B288
-------------- next part --------------
A non-text attachment was scrubbed...
Name: p_error.png
Type: image/png
Size: 16628 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-cs/attachments/20120508/c7a3a5a3/attachment-0004.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: p32.png
Type: image/png
Size: 16410 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-cs/attachments/20120508/c7a3a5a3/attachment-0005.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: p64.png
Type: image/png
Size: 20334 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-cs/attachments/20120508/c7a3a5a3/attachment-0006.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: p128.png
Type: image/png
Size: 27809 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-cs/attachments/20120508/c7a3a5a3/attachment-0007.png
More information about the CIG-CS
mailing list