[CIG-CS] Discontinuous viscosity

[Walter Landry]

Hello Everyone,

Back in May 2011, I ran a test for solving Stokes with discontinuous
viscosity using deal.II and Gamr, my finite difference code.

  http://www.geodynamics.org/pipermail/cig-cs/2011-May/000035.html

That highlighted a problem that all codes encounter with discontinuous
viscosities.  I think I now have a solution for this problem.  I have
implemented a variation of the Immersed Interface Method (IIM) in a
finite difference code.  The basic idea of IIM is to use normal finite
differences and then add correction terms arising from the jumps.

One interesting detail is that it is necessary to change the
fundamental variables, because I do not have expressions for the jumps
in velocity.  To that end, I use the augmented velocity, which is the
ordinary velocity multiplied by the viscosity.

In any case, to test the code I used SOLCX [1], a test used by other
groups to gauge how well a code handles variable viscosity [2] [3].
SOLCX is set in a unit box and has a sharp jump in viscosity at x=x_c.
For these tests, I used x_c=0.4, which is not aligned with any grid
boundary.  Moresi et al [3], with their finite element code, saw
errors in pressure of 80% with a 32x32 grid when the elements were not
aligned with the grid.  They do not report on convergence, though,
based on the work I did in May 2011, I suspect that it diverges.

In comparison, I am attaching a plot of the pressure solution and the
error at y=0.25 for a viscosity ratio of 10^10.  This method should
converge to first order on the boundary.  For this toy code I am using
Gauss-Seidel relaxation.  The number of iterations is largely
independent of the viscosity ratio.  For a 64x64 grid, using vx=vy=p=0
as the initial guess:

  viscosity ratio    # iterations
  1                  18415
  1e2                22685
  1e4                22764
  1e6                22765
  1e8                22765
  1e10               22765

As resolution improves, the number of iterations increases, but the
error in the pressure decreases.

  viscosity ratio=10^10
  h      # iterations  L_infinity(pressure)

  1/32   5769          .0184
  1/64   22765         .01562
  1/128  92946         .00246


As I mentioned, this is still a toy code.  It is hard coded for a
horizontal viscosity jump.  The solver is just Gauss-Seidel
relaxation, but it does use the same relaxation parameters as the full
parallel variable viscosity Stokes solver in Gamr.  It does end up
using a larger stencil than the straight Gamr code, but otherwise it
should apply to the full code without too much trouble.  If you are
interested, the code is available through mercurial with the command

  hg clone http://geodynamics.org/hg/cs/AMR/Discontinuous_Stokes/

One nice thing about this method is that it also works as a way of
embedding boundaries.  So I can use it to model a spherical earth
inside a rectilinear grid.

So the next step is to get this method working for more general
interfaces.  My plan is to embed SOLCX into a larger grid and then
rotate it.  After that, I will use an analytic solution for a circular
inclusion [4].  The last step is to put it into Gamr.

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] Preconditioned iterative methods for Stokes flow problems arising in
    computational geodynamics
    Dave A. May, Louis Moresi
    Physics of the Earth and Planetary Interiors 171 (2008) 33$(Q#|(B47

[3] The accuracy of finite element solutions of Stokes' flow with
    strongly varying viscosity
    Louis Moresi, Shijie Zhong, Michael Gurnis
    Physics of the Earth and Planetary Interiors 97 (1996) 83-94

[4] 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.png
Type: image/png
Size: 15116 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-cs/attachments/20120304/109b87da/attachment-0002.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: p_error.png
Type: image/png
Size: 16617 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-cs/attachments/20120304/109b87da/attachment-0003.png