Convergence issue scaling square root

1

Hi all,

I want to solve a diffusion equation with some scaling of the solution.

The scaling is theta = (c/S)^2

See example below:

from fenics import *

num_steps = 3  # number of time steps
dt = 0.0001  # time step size

# some scaling parameter
S = 2

# Create mesh and define function space
nx = 1000
mesh = IntervalMesh(nx, 0, 2)
V = FunctionSpace(mesh, 'CG', 1)

# Define initial value
c_0 = Expression('2 + x[0]', degree=0)
theta_0 = (c_0/S)**2  # theta = (c/S)**2
theta_n = project(theta_0, V)

# Define variational problem
theta = Function(V)
theta.assign(theta_n)  # initial guess
v = TestFunction(V)
f = Constant(0)

c = S*(theta + DOLFIN_EPS)**(1/2)  # c = S * sqrt(theta)
c_n = S*(theta_n + DOLFIN_EPS)**(1/2)  # c = S * sqrt(theta)

F = c*v*dx + dt*dot(grad(c), grad(v))*dx - (c_n + dt*f)*v*dx

bcs = [DirichletBC(V, Constant(0), "on_boundary")]

#####################################################
# Usual nonlinear problem formulation
du = TrialFunction(V)
a = derivative(F, theta, du)# L is the weak form
problem = NonlinearVariationalProblem(F, theta, bcs, J=a)
solver = NonlinearVariationalSolver(problem)

# Solver parameters
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-6
prm['newton_solver']['relative_tolerance'] = 1E-6
prm['newton_solver']['maximum_iterations'] = 20
#####################################################

t = 0
for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solver.solve()
    theta_n.assign(theta)

For some reason, there are some convergence issues:

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 4.124e+00 (tol = 1.000e-06) r (rel) = 1.000e+00 (tol = 1.000e-06)
  Newton iteration 1: r (abs) = 4.566e-02 (tol = 1.000e-06) r (rel) = 1.107e-02 (tol = 1.000e-06)
  Newton iteration 2: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 3: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 4: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 5: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 6: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 7: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 8: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 9: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 10: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 11: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 12: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 13: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 14: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 15: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 16: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 17: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 18: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 19: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
  Newton iteration 20: r (abs) = -nan (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)
Traceback (most recent call last):
  File "mwe_solution_to_power_solving.py", line 51, in <module>
    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.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

If we remove the line theta.assign(theta_n) # initial guess it works but the first iteration is very slow to converge.

Help is greatly appreciated!

Rem

2

I believe the source of the problem is that there is a jump discontinuity between the homogeneous Dirichlet data and the initial condition, which is producing oscillations in the numerical solution. These oscillations undershoot the zero value enforced at the boundary, resulting in a negative argument to the square root. If I apply homogeneous Neumann BCs (i.e., bcs = []) or use c_0 as the Dirichlet data, then there is no convergence problem. Is the given boundary condition what you want in your application, or just a placeholder for developing the code?

1 Like
3

Ideally I want users to be able to give a value for c_0, apply the scaled BC, solve for theta, and scale again to show users c.