Why does Solver not converge in specific cases?

1

Hello everyone,

I want to use dolfin_dg to solve the Euler equations. I have been working through the demos to check the convergence rate of the dg method in two dimensions which is done here. My problem is that I am trying to check the convergence with some other manufactured solutions and when using a different solution the solver does not obtain a satisfactory solution or does not converge. For example, when running the code with constant velocities and energy along the domain, the newton method arrives at a NaN for the residue in the first iteration. I attach a MWE and the output (where I erased all the iterations after the first NaN.)

Thanks in advance

from dolfin import *

from dolfin_dg import *

parameters["ghost_mode"] = "shared_facet"

parameters['form_compiler']['representation'] = 'uflacs'

parameters['form_compiler']["quadrature_degree"] = 8

ele_n = 32

p = 1

# Mesh and function space.

mesh = RectangleMesh(Point(0, 0), Point(.5*pi, .5*pi), ele_n, ele_n, 'right')

V = VectorFunctionSpace(mesh, 'DG', p, dim=4)

# Set up Dirichlet BC

gD = Expression(('sin(x[0])*sin(x[1]) + 2',

                'sin(x[0])*sin(x[1]) + 2',

                'sin(x[0])*sin(x[1]) + 2',

                'sin(x[0])*sin(x[1]) + 2'),

                element=V.ufl_element())

f = Expression(('sin(x[0])*cos(x[1]) + sin(x[1])*cos(x[0])',

                'sin(x[0])*cos(x[1]) + sin(x[1])*cos(x[0])',

                'sin(x[0])*cos(x[1]) + sin(x[1])*cos(x[0])',

                'sin(x[0])*cos(x[1]) + sin(x[1])*cos(x[0])'),

            element=V.ufl_element())

u, v = interpolate(gD, V), TestFunction(V)

bo = CompressibleEulerOperator(mesh, V, DGDirichletBC(ds, gD))

residual = bo.generate_fem_formulation(u, v) - inner(f, v)*dx

solve(residual == 0, u)

errorl2 = errornorm(gD, u, norm_type='l2', degree_rise=3)

print(errorl2)

The terminal returns:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F’.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 3.445e-04 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Traceback (most recent call last):
File “manufactured_sol_original.py”, line 36, in
solve(residual == 0, u)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 233, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 314, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

2

I actually noticed this yesterday myself.

I’m also not sure how this problem snuck into dolfin_dg's reproduction examples, since it’s not actually exhibited in the paper. Quite embarrassing!

I recommend looking at the Ringleb flow example as an extensive test of the solution techniques. You can find documentation in Section 4.2 here.

Nevertheless I’ll take a look at the Euler MMS example you’ve linked to and see if I can deduce the problem. Likely it’s not well defined in the (0, \pi / 2)^2 domain with influx BCs everywhere.