Hi Walter,<div>I'm very interested in what happen when the interface in not horizontal. Using some similar techniques, I've found very different results when the interface is slightly tilted. </div><div>Regards,</div>
<div>Sergio</div><div><br></div><div><br clear="all"><br><br><br>Sergio Zlotnik<br>LaCaN, Laboratori de Calcul Numeric<br>Dept. Matematica Aplicada III<br>UPC-Barcelona Tech<br>Jordi Girona 1-3 E-08034<br>Barcelona, Spain<br>
<br><a href="mailto:sergio.zlotnik@upc.edu" target="_blank">sergio.zlotnik@upc.edu</a><br>Ph: +34 93 401 7760<br>Fax: +34 93 401 1825<br>office: C2-203<br>
<br><br><div class="gmail_quote">2012/3/5 Walter Landry <span dir="ltr"><<a href="mailto:walter@geodynamics.org">walter@geodynamics.org</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hello Everyone,<br>
<br>
Back in May 2011, I ran a test for solving Stokes with discontinuous<br>
viscosity using deal.II and Gamr, my finite difference code.<br>
<br>
<a href="http://www.geodynamics.org/pipermail/cig-cs/2011-May/000035.html" target="_blank">http://www.geodynamics.org/pipermail/cig-cs/2011-May/000035.html</a><br>
<br>
That highlighted a problem that all codes encounter with discontinuous<br>
viscosities. I think I now have a solution for this problem. I have<br>
implemented a variation of the Immersed Interface Method (IIM) in a<br>
finite difference code. The basic idea of IIM is to use normal finite<br>
differences and then add correction terms arising from the jumps.<br>
<br>
One interesting detail is that it is necessary to change the<br>
fundamental variables, because I do not have expressions for the jumps<br>
in velocity. To that end, I use the augmented velocity, which is the<br>
ordinary velocity multiplied by the viscosity.<br>
<br>
In any case, to test the code I used SOLCX [1], a test used by other<br>
groups to gauge how well a code handles variable viscosity [2] [3].<br>
SOLCX is set in a unit box and has a sharp jump in viscosity at x=x_c.<br>
For these tests, I used x_c=0.4, which is not aligned with any grid<br>
boundary. Moresi et al [3], with their finite element code, saw<br>
errors in pressure of 80% with a 32x32 grid when the elements were not<br>
aligned with the grid. They do not report on convergence, though,<br>
based on the work I did in May 2011, I suspect that it diverges.<br>
<br>
In comparison, I am attaching a plot of the pressure solution and the<br>
error at y=0.25 for a viscosity ratio of 10^10. This method should<br>
converge to first order on the boundary. For this toy code I am using<br>
Gauss-Seidel relaxation. The number of iterations is largely<br>
independent of the viscosity ratio. For a 64x64 grid, using vx=vy=p=0<br>
as the initial guess:<br>
<br>
viscosity ratio # iterations<br>
1 18415<br>
1e2 22685<br>
1e4 22764<br>
1e6 22765<br>
1e8 22765<br>
1e10 22765<br>
<br>
As resolution improves, the number of iterations increases, but the<br>
error in the pressure decreases.<br>
<br>
viscosity ratio=10^10<br>
h # iterations L_infinity(pressure)<br>
<br>
1/32 5769 .0184<br>
1/64 22765 .01562<br>
1/128 92946 .00246<br>
<br>
<br>
As I mentioned, this is still a toy code. It is hard coded for a<br>
horizontal viscosity jump. The solver is just Gauss-Seidel<br>
relaxation, but it does use the same relaxation parameters as the full<br>
parallel variable viscosity Stokes solver in Gamr. It does end up<br>
using a larger stencil than the straight Gamr code, but otherwise it<br>
should apply to the full code without too much trouble. If you are<br>
interested, the code is available through mercurial with the command<br>
<br>
hg clone <a href="http://geodynamics.org/hg/cs/AMR/Discontinuous_Stokes/" target="_blank">http://geodynamics.org/hg/cs/AMR/Discontinuous_Stokes/</a><br>
<br>
One nice thing about this method is that it also works as a way of<br>
embedding boundaries. So I can use it to model a spherical earth<br>
inside a rectilinear grid.<br>
<br>
So the next step is to get this method working for more general<br>
interfaces. My plan is to embed SOLCX into a larger grid and then<br>
rotate it. After that, I will use an analytic solution for a circular<br>
inclusion [4]. The last step is to put it into Gamr.<br>
<br>
Cheers,<br>
Walter Landry<br>
<br>
<br>
[1] Analytic solutions for Stokes’ flow with lateral variations in viscosity<br>
Shijie Zhong<br>
Geophys. J. Int.(1996) 124, 18-28<br>
Section 2.3<br>
<br>
[2] Preconditioned iterative methods for Stokes flow problems arising in<br>
computational geodynamics<br>
Dave A. May, Louis Moresi<br>
Physics of the Earth and Planetary Interiors 171 (2008) 33�#|47<br>
<br>
[3] The accuracy of finite element solutions of Stokes' flow with<br>
strongly varying viscosity<br>
Louis Moresi, Shijie Zhong, Michael Gurnis<br>
Physics of the Earth and Planetary Interiors <a href="tel:97%20%281996%29%2083-94" value="+19719968394">97 (1996) 83-94</a><br>
<br>
[4] Analytical solutions for deformable elliptical inclusions<br>
in general shear<br>
Daniel W. Schmid and Yuri Yu. Podladchikov<br>
Geophys. J. Int. (2003) 155, 269�#|288<br>
<br>_______________________________________________<br>
CIG-CS mailing list<br>
<a href="mailto:CIG-CS@geodynamics.org">CIG-CS@geodynamics.org</a><br>
<a href="http://geodynamics.org/cgi-bin/mailman/listinfo/cig-cs" target="_blank">http://geodynamics.org/cgi-bin/mailman/listinfo/cig-cs</a><br>
<br></blockquote></div><br></div>