### Newton solver diverges

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?

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?

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.

https://arxiv.org/abs/1609.02190

Hope that helps.

### 1 Answer

```
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.

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.

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
```

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.