### Newton solver diverges

260
views
0
7 months ago by
HI all

I'm struggling with mesh convergence of the problem described here: https://www.allanswered.com/post/xnnob/bc-for-a-hyperbolic-1st-order-system-of-pde/
I'm using the solution proposed by pf4d.

When I increase the number of elements in the mesh I would expect the solutions to converge. However, the min and max values fluctuate widely.
I have solved the problem for increasing mesh densitiy N and printed the min and max of the solution. This is the output:
N=  32
min =  -14.213930736628392  max =  9.046530167293112
N=  64
min =  -2.2292666922532374  max =  2.5227212758820046
N=  128
min =  -499.57598020660737  max =  460.7371101576128
N=  256
min =  -6.909612880745566  max =  6.33138530677423
N=  512
min =  -5.53777342589658  max =  3.98577741119527
N=  1024
min =  -30.5893553449149  max =  23.895255532377046
N=  2048
....
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.​

If I use econd order elements instead of first order elements,
P1      = FiniteElement('CG', triangle,2)​
the solutions clearly diverge:
N=  32
min =  -2.5901446402149064  max =  1.6248883001122068
N=  64
min =  -61.72313541884393  max =  98.2433484043985
N=  128
min =  -133.96435807556034  max =  176.622721145695
N=  256
min =  -229.42673119333665  max =  238.81759600560238
N=  512
min =  -992.7646840741731  max =  968.4059946564718
N=  1024
...
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.​
I guess that the error message comes because the mesh is too fine and /or numbers are getting too big.

Should I use another solver? Which one might give a better convergence? Any other proposals what I could do?
Community: FEniCS Project
Hyperbolic equations are unstable using standard finite-element function spaces, and require augmentation in order to be stated as a well-posed discrete system.

The system of equations you wish to solve can be re-written with velocity vector $\underline{u} = [u\ v]^{\intercal}$ as
\begin{align} \nabla \cdot \underline{u} &= 0 \\ \frac{\partial u}{\partial y} &= \frac{\partial v}{\partial x} \end{align},
The first equation is conservation of mass for an incompressible fluid while the second enforces symmetry of the velocity-gradient tensor
$\nabla \underline{u} = \begin{bmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{bmatrix}.$
I am curious to know what physical system you are attempting to approximate.

Anyway, I suspect that you are going to need to stabilize the numerical system of equations.  I believe the reason you get reduced-magnitude oscillations with lower mesh resolution is due to the numerical diffusion coarser meshes provide.
written 7 months ago by pf4d
The funny thing is that Mathematica solves this problem as it is without any problems (except that it is slow and tends to crash)...

I'm trying to solve a simplified form of the stress balance at a calving front. In general the stress balance is
$\nabla\cdot \sigma = f$
with the stress tensor $\sigma$ and the force vector f. I neglect one of the horizontal directions and introduce a deviatoric stress tensor, such that $\sigma = S + PI$. Including non-compressibility ($S_{xx} = -S_{zz}$) this gives
$\sigma = \begin{pmatrix} S_{xx} + \rho g z & S_{xz} \\ S_{xz} & - S_{xx} + \rho g z \end{pmatrix}$
Then I define $A=S_{xx}$ and $B=S_{xz}$ and get the equations
$\partial_x A + \partial_z B = 0$
$\partial_x B - \partial_z A = 0$
I tried to do something similar to what you described, but because of the incompressibility I have this minus sign in the second equation.

How do I stabilize the system? By rewriting it as a relaxation problem?
written 7 months ago by Ta S
By specifying the pressure as hydrostatic, you lose one unknown and must therefore incorporate incompressibility into the momentum system explicitly.  I have not seen this done before, that I can recall.

Why don't you just solve the plane-strain Stokes equations?  It is quite efficient, with a plethora of documentation,

https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/stokes-stabilized/python/documentation.html

for example.
written 7 months ago by pf4d
All examples I saw of the Stokes equation solve for a vector u of displacement, that means I need a constitutve relation to convert stress to strain, i.e. linear elasticity or sth more complicated. Since this is not trivial for ice I want to avoid it.
written 7 months ago by Ta S
Glen's flow law is not difficult to implement, see e.g.

https://arxiv.org/abs/1609.02190

Hope that helps.
written 7 months ago by pf4d
I know that, what I mean is that the choice of rheology is not trivial. Which exponent should I use in Glen's flow law?  n=3 is standard but not necessarily the best choice for my application. It's an additional assumption I need to make and one that I don't actually need at the moment because I'm only interested in stresses and not in flow.
written 7 months ago by Ta S
Regardless, I think solving for the deviatoric stress distribution like this is a neat idea.
written 7 months ago by pf4d

0
7 months ago by
I have restated the problem here and have applied some artificial diffusion in the direction of stress : see if the output from this script makes sense to you:

from fenics import *
from mshr   import generate_mesh, Rectangle

# Define domain and create mesh
H       = 0.1
L       = 1.0
domain  = Rectangle(Point(0, 0), Point(L, H))
mesh    = generate_mesh(domain, 400)

# Define the function space
P1      = FiniteElement('CG', triangle, 1)
element = MixedElement([P1, P1])
V       = FunctionSpace(mesh, element)

# Define test and trial functions
psi    = TestFunction(V)
phi    = TrialFunction(V)
u      = Function(V)
h      = CellSize(mesh)
f      = Constant((0,0))

# Define boundaries
def surface(x, on_boundary):  return on_boundary and near(x[1], 0)
def glacier(x, on_boundary):  return on_boundary and near(x[0], L)
def front(x, on_boundary):    return on_boundary and near(x[0], 0)

# Define front stress
front_normal_stress = Expression('-x[1]', degree=1)

# Define boundary conditions
bc_surface = DirichletBC(V, Constant((0, 0)), surface)
bc_glacier = DirichletBC(V, Constant((0, 0)), glacier)
bc_front   = DirichletBC(V.sub(0), front_normal_stress, front)
bc         = [bc_glacier, bc_surface, bc_front]

# SUPG intrinsic time parameter :
tau   = 10*h

# the deviatoric stress tensor :
def sigma(u):
return as_matrix([[u[0],  u[1]],
[u[1], -u[0]]])

# the differential operator :
def Lu(u): return div(sigma(u))

# Define the variational problem :
F = inner(Lu(u), psi)*dx - inner(f, psi)*dx \
+ inner(Lu(psi), tau*Lu(u)) * dx \
- inner(Lu(psi), tau*f) * dx

# define the Jacobian :
J = derivative(F, u, phi)

# compute solution using Newton's method :
solve(F==0, u, bc, J=J)

print u.vector().min(), u.vector().max()

File('u.pvd') << u
​

Note that I have adjusted "tau" to be dependant on the mesh size, and found by experimentation that a factor of 10 produces a smooth result ...  You should experiment with this value.  For mathematical background regarding the technique, look up "Galerkin/least-squares" stabilization, or read the paper:

Thomas J. R. HUGHES, Leopoldo E. FRANCA and Gregory M. HULBERT (1988): A NEW FINITE ELEMENT FORMULATION FOR COMPUTATIONAL FLUID DYNAMICS: VIII. GALERKIN/LEAST-SQUARES METHOD FOR ADVECTIVE-DIFFUSIVE EQUATIONS.

Note:
Usually, water pressure is applied over the terminus, but as you have specified hydrostatic internal pressure and eliminated the forcing term $f$ from the momentum balance, you are specifying deviations from the mean pressure.

Hm... it looks similar to what I would expect, but the stress maxima should be close to the front and not at the glacier boundary - that makes no sense.
Thanks for all your effort. I've decided to go back to Mathematica, because it solves the equations easily even though working with the solutions is a bit of a hassle.
written 7 months ago by Ta S
Yes, this is because the boundary conditions are wrong; I was never certain what you were trying to do with the BCS.  The above code specifies the surface at $x=0$ -- I assume you want this to be the bed.

Below, I've imposed zero-deviatoric stress conditions everywhere except on the "front" where I have imposed a linearly-increasing with depth essential boundary on the $x$ component of deviatoric stress.  The $y$ component of deviatoric stress along the "front" boundary is allowed to vary freely such that the PDE is satisfied, i.e., so that $\nabla \cdot \sigma = \underline{0}$, where $\sigma$ is the rank-two deviatoric-stress tensor under hydrostatic-pressure and plane-strain approximations.

from fenics import *
from mshr   import generate_mesh, Rectangle

# Define domain and create mesh
H       = 0.1
L       = 1.0
domain  = Rectangle(Point(0, -H), Point(L, 0))
mesh    = generate_mesh(domain, 400)

# Define the function space
P1      = FiniteElement('CG', triangle, 1)
element = MixedElement([P1, P1])
V       = FunctionSpace(mesh, element)

# Define test and trial functions
psi    = TestFunction(V)
phi    = TrialFunction(V)
u      = Function(V)
h      = CellSize(mesh)
f      = Constant((0,0))

# Define boundaries
def surface(x, on_boundary):  return on_boundary and near(x[1], 0)
def glacier(x, on_boundary):  return on_boundary and near(x[0], L)
def front(x, on_boundary):    return on_boundary and near(x[0], 0)
def bed(x, on_boundary):      return on_boundary and near(x[1], -H)

# Define front stress
front_normal_stress = Expression('-x[1]', degree=1)

# Define boundary conditions
bc_surface = DirichletBC(V, Constant((0, 0)), surface)
bc_glacier = DirichletBC(V, Constant((0, 0)), glacier)
bc_front   = DirichletBC(V.sub(0), front_normal_stress, front)
bc_bed     = DirichletBC(V, Constant((0, 0)), bed)
bc         = [bc_glacier, bc_surface, bc_front, bc_bed]

# SUPG intrinsic time parameter :
tau   = 10*h

# the deviatoric-stress tensor :
def sigma(u):
return as_matrix([[u[0],  u[1]],
[u[1], -u[0]]])

# the differential operator :
def Lu(u): return div(sigma(u))

# Define the variational problem :
F = inner(Lu(u), psi)*dx - inner(f, psi)*dx \
+ inner(Lu(psi), tau*Lu(u)) * dx \
- inner(Lu(psi), tau*f) * dx

# define the Jacobian :
J = derivative(F, u, phi)

# compute solution using Newton's method :
solve(F==0, u, bc, J=J)

print u.vector().min(), u.vector().max()

File('u.pvd') << u
​
written 7 months ago by pf4d