[CIG-CS] Variable viscosity Stokes solver

[Walter Landry]

Moving this discussing to the cig-cs list.

Wolfgang Bangerth <bangerth at math.tamu.edu> wrote:
> 
>> Taras Gerya manages to get up to 10^6 viscosity jumps for this problem
>> by using continuations.  However, when doing geophysical runs, he
>> never needs it.
>> 
>> Rhea has some large viscosity variations in their simulations, but I
>> do not think that their element-to-element variation is this large.
> 
> I don't think so either. In Rhea, the viscosity is determined by the 
> temperature and flow field, and a variation in eta so large would require an 
> element-to-element variation in temperature or strain rate that is also very 
> large. I do not believe that you would be able to get such variations in 
> solution fields of this magnitude from most reasonable finite element methods 
> -- you need to stabilize methods to resolve variations this large, and this 
> tends to spread the variation out over several cells, which would then also 
> spread out the variation in viscosity. 
> 
> In other words, requiring a solver to be able to deal with such large 
> viscosity variations is not something you'd likely encounter in a code in 
> which the viscosity is not a parameter but a function of other variables.

In problems with faulting, we can get arbitrarily small regions with
arbitrarily large viscosity jumps.  We can apply a limiter in order to
get a convergent algorithm, but the viscosity gradient will still be
very sharp.  For example, in Gale's dike example, the viscosity
variation is (I think) 10^5 element-to-element.  That is why I ended
up using a direct solver for it.

>> Dave May's solver handles 10^6, but he aligns the boundary of the
>> falling box with the edges of the element.
> 
> That actually leads to another interesting question: if the interface is not 
> aligned, what do you use to integrate the local matrices and vectors? If your 
> coefficient, Gauss formula are no longer useful. Presumably an iterated Gauss 
> formula, or something else of low order would be more appropriate.

Gale partitions the element up with Voronoi cells centered around each
particle.  The viscosity is assumed smooth on those cells.  So the
integration is essentially done by using the particles as integration
points.

I have also seen good results for other problems by using a higher
order field (e.g. fourth order) for material properties.  I do not
know if that will work for these problems.

Cheers,
Walter Landry
walter at geodynamics.org