[aspect-devel] Free-Surface in 3d

[Wolfgang Bangerth]

On 07/22/2015 12:34 PM, Ian Rose wrote:
>
> If we consider the typical size of entries in the normal system matrix, we get
> something like \eta L, where \eta is the viscosity and L is a lengthscale for
> the cells.  For free surface problems, the surface should be stabilized.  The
> stabilization term is typically of order \rho g \Delta t L^2.  The ratio of
> the size of the stabilization term to the normal matrix elements is then \rho
> g L \Delta t / \eta. For most problems of geophysical relevance this term will
> be on the small (but not negligible, since we need it for stabilization!)
> side, something like 10^{-1} - 10^{-2},
>
> However, the convection-box-3d parameter file has some odd values which change
> the story.  In particular, the value for gravity is large (10^{16}!!) and the
> value for viscosity is small (1).  So the stabilization term absolutely
> dominates the matrix entries by around nine or ten orders of magnitude.  You
> can imagine that that would be bad.  I would note that the unusual values for
> the convection-box parameter files have caused problems in the past, see, for
> example issue #94.  Perhaps Wolfgang or Timo would have some further insight here.

I haven't followed the development of the free surface terms from a 
mathematical perspective (I only looked at code issues), so I have nothing of 
real substance to offer other than the following question:

Stabilization terms are usually added to existing terms in the equation. As 
usual, the terms you add need to have the same physical units as the existing 
terms. Is this not the case here? Or is it the case that you have a 
stabilization parameter that has physical dimensions? One usually avoids such 
scaling problems by defining the stabilization term as the product of
- a dimensionless number
- a product of physical constants (such as the viscosity, density, etc) so
   that the overall term has the correct units consistent with the other
   terms in the equation
- the stabilization itself (e.g., the Laplacian of some quantity)

This way, one only has to tune the dimensionless first factor, and the whole 
thing will continue to work in other settings.

Best
  Wolfgang

-- 
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
                                 www: http://www.math.tamu.edu/~bangerth/