# Modified Euler Equations with dolfin_dg

Hello everyone, I am using dolfin_dg to perform some simulations of modified euler equations. basically what I want to solve is the euler equations with a modified pressure and a polytropic relation substituting the energy equation.

\frac{d \vec q_j}{dt} + \vec\nabla_i \textit{F}_{i,j} = 0

where

q_j = (n,n m u_1, n m u_2)^T\\ F_{i,j}=(n\vec u, m n \vec( u \otimes\vec u)+ T n I_{2 \times 2})^T

with T a constant.

I tried to supply the flux as daughter class inheriting from the HyperbolicOperator class:

class PlasmaOperator(HyperbolicOperator):

def __init__(self, mesh, V, bcs):

dim = mesh.geometric_dimension()

def F_c(U):

n_i = U

n_m_u = as_vector([U[j] for j in range(1, dim+1)])

u = n_m_u/(m_i*n_i)

p = n_i*electron_temp

inertia = n_i*ufl.outer(u, u) + p*Identity(dim)

res = ufl.as_tensor([n_i*u,

*[inertia[d, :] for d in range(dim)]])
return res

def alpha(U, n):

n_i = U

n_m_u = as_vector([U[j] for j in range(1, dim+1)])

u = n_m_u/(m_i*n_i+DOLFIN_EPS)

c = sqrt(T/m_i)

lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]

return lambdas

HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, LocalLaxFriedrichs(alpha))


The boundary conditions I am trying to implement are basically an Dirichlet in the left side of the box and free outflow anywhere, while the initial condition should be zero density and velocity everywhere, this generates problems as you get 0/0, anyway the zero density is not a strong constraint as I could change it. Here’s how I implement them:

gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))

class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x < DOLFIN_EPS

class Outlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x > DOLFIN_EPS

bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)

INLET, OUTFLOW = 1, 2

inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)
outlet = Outlet()
outlet.mark(bound_mshfunc, OUTFLOW)

ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)

bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGDirichletBC(ds(OUTFLOW),gD2)]


However when I solve the problem the residual from Newton’s solver increases steadily and diverges.

Here is my full code:

#%%
import matplotlib.pyplot as plt
import ufl
from dolfin import *
from dolfin_dg import *

parameters["ghost_mode"] = "shared_facet"

# Mesh and function space.

N_ele = 10
p = 2
mesh = UnitSquareMesh(N_ele,N_ele,'left/right')
V = VectorFunctionSpace(mesh, 'DG', p, dim=3)
du = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)

#Constants
electron_temp = 20 #T
m_i = 1
e = 1.6

n0 = Constant(1)
u0 = Constant(1)
v0 = Constant(0.0)

gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))
gD2 = as_vector((0.1, 0, 0))

u_vec_0 = Constant((n0, m_i*n0*u0, m_i*n0*v0))
u_vec = interpolate(u_vec_0,V)

T = 0.2
h = 1/N_ele
CFL = 0.6
dt_n = CFL*h/(p+1)**2
N = int(T/dt_n)
dt = Constant(dt_n)
un = Function(V)

###################################################

class PlasmaOperator(HyperbolicOperator):

def __init__(self, mesh, V, bcs):

dim = mesh.geometric_dimension()

def F_c(U):
n_i = U
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i)
p = n_i*electron_temp

inertia = n_i*ufl.outer(u, u) + p*Identity(dim)
res = ufl.as_tensor([n_i*u,
*[inertia[d, :] for d in range(dim)]])

return res

def alpha(U, n):
n_i = U
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i+DOLFIN_EPS)
c = sqrt(electron_temp/m_i)
lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]
return lambdas

HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, LocalLaxFriedrichs(alpha))

###################################################

class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x < DOLFIN_EPS

class Outlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x > DOLFIN_EPS

bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)

INLET, OUTFLOW = 1, 2

inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)
outlet = Outlet()
outlet.mark(bound_mshfunc, OUTFLOW)

ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)

bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGDirichletBC(ds(OUTFLOW),gD2)]

po = PlasmaOperator(mesh, V, bcs= bcs)
residual = po.generate_fem_formulation(u_vec, v)
residual += dot(u_vec - un, v)*dx

J = derivative(residual, u_vec, du)

class Problem(NonlinearProblem):
def F(self, b, x):
assemble(residual, tensor=b)

def J(self, A, x):
assemble(J, tensor=A)

solver = NewtonSolver()

for j in range(N):
un.assign(u_vec)
solver.solve(Problem(), u_vec.vector())



1 Like

This is a tricky one. I’ll have a think about it and get back to you.

import matplotlib.pyplot as plt
import ufl
from dolfin import *
from dolfin_dg import *

parameters["ghost_mode"] = "shared_facet"

# Mesh and function space.

N_ele = 10
p = 2
mesh = UnitSquareMesh(N_ele,N_ele,'left/right')
V = VectorFunctionSpace(mesh, 'DG', p, dim=3)
du = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)

#Constants
electron_temp = 20 #T
m_i = 1
e = 1.6

n0 = Constant(1)
u0 = Constant(1)
v0 = Constant(0.0)

# These BCs are not continuous... be careful
gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))
gD2 = as_vector((0.1, 0, 0))

u_vec_0 = Constant((n0, m_i*n0*u0, m_i*n0*v0))
u_vec = interpolate(u_vec_0,V)

T = 0.2
h = 1/N_ele
CFL = 0.6
dt_n = CFL*h/(p+1)**2
N = int(T/dt_n)
dt = Constant(dt_n)
un = Function(V)

###################################################

class PlasmaOperator(HyperbolicOperator):

def __init__(self, mesh, V, bcs):

dim = mesh.geometric_dimension()

def F_c(U):
n_i = U
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i)
p = n_i*electron_temp

inertia = n_i*ufl.outer(u, u) + p*Identity(dim)
res = ufl.as_tensor([n_i*u,
*[inertia[d, :] for d in range(dim)]])

return res

def alpha(U, n):
n_i = U
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i+DOLFIN_EPS)
c = sqrt(electron_temp/m_i)
lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]
return lambdas

HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, LocalLaxFriedrichs(alpha))

###################################################

class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0.0)

bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)

INLET, OUTFLOW = 1, 2

# Mark the whole exterior as outflow
AutoSubDomain(lambda x, on_boundary: on_boundary).mark(bound_mshfunc, OUTFLOW)

# Then overwrite the inflow markers
inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)

ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)

bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGDirichletBC(ds(OUTFLOW),gD2)]

po = PlasmaOperator(mesh, V, bcs= bcs)

# Don't forget dt in the temporal discretisation
residual = dt*po.generate_fem_formulation(u_vec, v)
residual += dot(u_vec - un, v)*dx

J = derivative(residual, u_vec, du)

class Problem(NonlinearProblem):
def F(self, b, x):
assemble(residual, tensor=b)

def J(self, A, x):
assemble(J, tensor=A)

solver = NewtonSolver()

for j in range(N):
un.assign(u_vec)
solver.solve(Problem(), u_vec.vector())
XDMFFile("u.xdmf").write_checkpoint(u_vec.sub(0), "u", append=j!=0, time_step=float(j))


In your previous attempt you didn’t capture all of the exterior boundaries. Note the top left and bottom left. I fixed this with:

# Mark the whole exterior as outflow
AutoSubDomain(lambda x, on_boundary: on_boundary).mark(bound_mshfunc, OUTFLOW)
# Then overwrite the inflow markers
inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)


You forgot dt in your temporal discretisation

# Don't forget dt in the temporal discretisation
residual = dt*po.generate_fem_formulation(u_vec, v)


Be very careful with your BCs. They are not continuous

# These BCs are not continuous... be careful
gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))
gD2 = as_vector((0.1, 0, 0))


I find it useful to use a finite volume method (i.e. p=0) initially to debug instability issues. FVMs are typically very diffusive and robust and at least let you see a “solution”. I used this to see the missing exterior facets problem.

Here is a snapshot of what I see of the “corrected” (but probably still not physically correct) n_h solution. Looks like the ilposed BCs and initial state is propagating error through time. 1 Like

Hi @nate, thanks for the fast reply.

That was very helpful, I had not realized that the boundaries were not correctly marked. On the other hand it is true that the problem is not well posed, I’ll implement continuous boundary conditions and be more careful with the initial conditions.

Thanks a lot for the tips.

1 Like

Hi @nate ,

I have made some changes in the boundary conditions, basically filling an empty tube with an inflow boundary condition in the left wall and free outflow everywhere else.

gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))
gD2 = as_vector((0, 0, 0))

class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0.0)

bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)

INLET, OUTFLOW = 1, 2

# Mark the whole exterior as outflow
AutoSubDomain(lambda x, on_boundary: on_boundary).mark(bound_mshfunc, OUTFLOW)

# Then overwrite the inflow and outflow markers
inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)

ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)

bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGNeumannBC(ds(OUTFLOW), gD2)]



The solution looks pretty much as I expect until some point in time when the solver does not converge. In fact the solver gets stuck at a relative tolerance of \sim 1.6e-6. Might there be a reason why a higher accuracy cannot be reached?

Other concerns I have about the behaviour of the solver are:

• In all timeteps in the first iteration of the solver the residual increases before decreasing and converging to the solution I don’t know if this is common or if I should check if my residual is well defined.
• In the first timestep the solver takes much longer to converge than in the other ones (18 iterations vs 3 or 4 in the subsequent ones).

Although this two last points were surprising the convergence rate does not look too bad (I think). I attach below the meaningfull part of the code.

import matplotlib.pyplot as plt
import ufl
from dolfin import *
from dolfin_dg import *

parameters["ghost_mode"] = "shared_facet"

# Mesh and function space.
x1, x2 = 0, 3
y1, y2 = -0.5, 0.5
N_ele = 10
p = 1
mesh = RectangleMesh(Point(x1,y1), Point(x2,y2),int(x2-x1)*N_ele,int(y2-y1)*N_ele,'left/right')
V = VectorFunctionSpace(mesh, 'DG', p, dim=3)
du = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)

electron_temp = 20
m_i = 4.2e-7
e = 1.6e-19

n0 = Constant(7e18)
u0 = Constant(1.4e5)
v0 = Constant(0)

gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))
gD2 = as_vector((0, 0, 0))

u_vec_0 = Expression(('0.00001*n0','0.00001* m_i*n0*u0', '0.00001*m_i*n0*v0'), n0=n0, m_i=m_i, u0=u0,v0=v0, element= V.ufl_element())
u_vec = interpolate(u_vec_0,V)

# time step
T = 5e-5
h = 1/N_ele
dt_n = 0.02*T
N = int(T/dt_n)
dt = Constant(dt_n)
un = Function(V)

###################################################
#Define Boundaries and boundary conditions
class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0.0)

bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)

INLET, OUTFLOW = 1, 2

# Mark the whole exterior as outflow
AutoSubDomain(lambda x, on_boundary: on_boundary).mark(bound_mshfunc, OUTFLOW)

# Then overwrite the inflow and outflow markers
inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)

ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)

bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGNeumannBC(ds(OUTFLOW), gD2)]

#############################################################

# Generate the DGFEM formulation for the operator

class Operator(HyperbolicOperator):

def __init__(self, mesh, V, bcs):

dim = mesh.geometric_dimension()

def F_c(U):
n_i = U
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i)
p = n_i*electron_temp

inertia = n_i*ufl.outer(u, u) + p*Identity(dim)
res = ufl.as_tensor([n_i*u,
*[inertia[d, :] for d in range(dim)]])

return res

def alpha(U, n):
n_i = U
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i)
c = sqrt(electron_temp/m_i)
lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]
return lambdas

HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, LocalLaxFriedrichs(alpha))

###################################################

po = Operator(mesh, V, bcs= bcs)
residual = dt*po.generate_fem_formulation(u_vec, v)
residual += dot(u_vec - un, v)*dx

f = as_vector((0,0,0))

residual += - dot(f,v)*dx(mesh)

J = derivative(residual, u_vec, du)

class Problem(NonlinearProblem):
def F(self, b, x):
assemble(residual, tensor=b)

def J(self, A, x):
assemble(J, tensor=A)

solver = NewtonSolver()
solver.parameters["relative_tolerance"] = 1e-6

# Run the simulation until elapsed time = T

for j in range(N):
un.assign(u_vec)
print(" Step: ", j+1, "// t =", (j+1)*dt_n)
solver.solve(Problem(), u_vec.vector())

XDMFFile("u.xdmf").write_checkpoint(u_vec.sub(0), "u", append=j!=0, time_step=float(j))



And the output:

 Step:  1 // t = 2.0000000000000003e-06
Newton iteration 0: r (abs) = 3.068e+22 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 3.031e+26 (tol = 1.000e-10) r (rel) = 9.879e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.260e+26 (tol = 1.000e-10) r (rel) = 4.106e+03 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 5.371e+25 (tol = 1.000e-10) r (rel) = 1.751e+03 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 2.344e+25 (tol = 1.000e-10) r (rel) = 7.642e+02 (tol = 1.000e-06)
Newton iteration 5: r (abs) = 1.040e+25 (tol = 1.000e-10) r (rel) = 3.391e+02 (tol = 1.000e-06)
Newton iteration 6: r (abs) = 4.679e+24 (tol = 1.000e-10) r (rel) = 1.525e+02 (tol = 1.000e-06)
Newton iteration 7: r (abs) = 2.127e+24 (tol = 1.000e-10) r (rel) = 6.933e+01 (tol = 1.000e-06)
Newton iteration 8: r (abs) = 9.736e+23 (tol = 1.000e-10) r (rel) = 3.174e+01 (tol = 1.000e-06)
Newton iteration 9: r (abs) = 4.473e+23 (tol = 1.000e-10) r (rel) = 1.458e+01 (tol = 1.000e-06)
Newton iteration 10: r (abs) = 2.051e+23 (tol = 1.000e-10) r (rel) = 6.686e+00 (tol = 1.000e-06)
Newton iteration 11: r (abs) = 9.301e+22 (tol = 1.000e-10) r (rel) = 3.032e+00 (tol = 1.000e-06)
Newton iteration 12: r (abs) = 4.083e+22 (tol = 1.000e-10) r (rel) = 1.331e+00 (tol = 1.000e-06)
Newton iteration 13: r (abs) = 1.651e+22 (tol = 1.000e-10) r (rel) = 5.380e-01 (tol = 1.000e-06)
Newton iteration 14: r (abs) = 5.436e+21 (tol = 1.000e-10) r (rel) = 1.772e-01 (tol = 1.000e-06)
Newton iteration 15: r (abs) = 1.016e+21 (tol = 1.000e-10) r (rel) = 3.313e-02 (tol = 1.000e-06)
Newton iteration 16: r (abs) = 4.279e+19 (tol = 1.000e-10) r (rel) = 1.395e-03 (tol = 1.000e-06)
Newton iteration 17: r (abs) = 7.052e+16 (tol = 1.000e-10) r (rel) = 2.299e-06 (tol = 1.000e-06)
Newton iteration 18: r (abs) = 3.782e+11 (tol = 1.000e-10) r (rel) = 1.233e-11 (tol = 1.000e-06)
Newton solver finished in 18 iterations and 18 linear solver iterations.
Step:  2 // t = 4.000000000000001e-06
Newton iteration 0: r (abs) = 8.857e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 3.296e+21 (tol = 1.000e-10) r (rel) = 3.721e+04 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.716e+20 (tol = 1.000e-10) r (rel) = 1.937e+03 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 7.165e+17 (tol = 1.000e-10) r (rel) = 8.089e+00 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.392e+13 (tol = 1.000e-10) r (rel) = 1.572e-04 (tol = 1.000e-06)
Newton iteration 5: r (abs) = 1.606e+08 (tol = 1.000e-10) r (rel) = 1.813e-09 (tol = 1.000e-06)
Newton solver finished in 5 iterations and 5 linear solver iterations.
Step:  3 // t = 6.000000000000001e-06
Newton iteration 0: r (abs) = 5.862e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.062e+21 (tol = 1.000e-10) r (rel) = 1.812e+04 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.258e+19 (tol = 1.000e-10) r (rel) = 3.853e+02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.313e+16 (tol = 1.000e-10) r (rel) = 2.240e-01 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 5.207e+09 (tol = 1.000e-10) r (rel) = 8.884e-08 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  4 // t = 8.000000000000001e-06
Newton iteration 0: r (abs) = 4.886e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 5.217e+20 (tol = 1.000e-10) r (rel) = 1.068e+04 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 5.959e+18 (tol = 1.000e-10) r (rel) = 1.220e+02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 9.106e+14 (tol = 1.000e-10) r (rel) = 1.864e-02 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.738e+08 (tol = 1.000e-10) r (rel) = 3.558e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  5 // t = 1.0000000000000003e-05
Newton iteration 0: r (abs) = 4.303e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 3.077e+20 (tol = 1.000e-10) r (rel) = 7.151e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.205e+18 (tol = 1.000e-10) r (rel) = 5.125e+01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.240e+14 (tol = 1.000e-10) r (rel) = 2.883e-03 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.668e+08 (tol = 1.000e-10) r (rel) = 3.876e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  6 // t = 1.2000000000000002e-05
Newton iteration 0: r (abs) = 3.891e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.019e+20 (tol = 1.000e-10) r (rel) = 5.187e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 9.986e+17 (tol = 1.000e-10) r (rel) = 2.566e+01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 2.550e+13 (tol = 1.000e-10) r (rel) = 6.553e-04 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.582e+08 (tol = 1.000e-10) r (rel) = 4.066e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  7 // t = 1.4000000000000001e-05
Newton iteration 0: r (abs) = 3.573e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.420e+20 (tol = 1.000e-10) r (rel) = 3.975e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 5.166e+17 (tol = 1.000e-10) r (rel) = 1.446e+01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 6.915e+12 (tol = 1.000e-10) r (rel) = 1.935e-04 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.618e+08 (tol = 1.000e-10) r (rel) = 4.529e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  8 // t = 1.6000000000000003e-05
Newton iteration 0: r (abs) = 3.311e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.050e+20 (tol = 1.000e-10) r (rel) = 3.170e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.936e+17 (tol = 1.000e-10) r (rel) = 8.867e+00 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 2.284e+12 (tol = 1.000e-10) r (rel) = 6.897e-05 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.684e+08 (tol = 1.000e-10) r (rel) = 5.086e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  9 // t = 1.8000000000000004e-05
Newton iteration 0: r (abs) = 3.085e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 8.050e+19 (tol = 1.000e-10) r (rel) = 2.609e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.790e+17 (tol = 1.000e-10) r (rel) = 5.800e+00 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 8.722e+11 (tol = 1.000e-10) r (rel) = 2.827e-05 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.596e+08 (tol = 1.000e-10) r (rel) = 5.172e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  10 // t = 2.0000000000000005e-05
Newton iteration 0: r (abs) = 2.884e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 6.349e+19 (tol = 1.000e-10) r (rel) = 2.202e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.151e+17 (tol = 1.000e-10) r (rel) = 3.991e+00 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 3.718e+11 (tol = 1.000e-10) r (rel) = 1.289e-05 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.575e+08 (tol = 1.000e-10) r (rel) = 5.460e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  11 // t = 2.2000000000000003e-05
Newton iteration 0: r (abs) = 2.699e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 5.121e+19 (tol = 1.000e-10) r (rel) = 1.898e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 7.718e+16 (tol = 1.000e-10) r (rel) = 2.860e+00 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.725e+11 (tol = 1.000e-10) r (rel) = 6.392e-06 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.578e+08 (tol = 1.000e-10) r (rel) = 5.846e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  12 // t = 2.4000000000000004e-05
Newton iteration 0: r (abs) = 2.525e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 4.204e+19 (tol = 1.000e-10) r (rel) = 1.665e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 5.350e+16 (tol = 1.000e-10) r (rel) = 2.119e+00 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 8.551e+10 (tol = 1.000e-10) r (rel) = 3.387e-06 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.567e+08 (tol = 1.000e-10) r (rel) = 6.206e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  13 // t = 2.6000000000000005e-05
Newton iteration 0: r (abs) = 2.359e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 3.500e+19 (tol = 1.000e-10) r (rel) = 1.483e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 3.807e+16 (tol = 1.000e-10) r (rel) = 1.614e+00 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 4.462e+10 (tol = 1.000e-10) r (rel) = 1.891e-06 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.558e+08 (tol = 1.000e-10) r (rel) = 6.604e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  14 // t = 2.8000000000000003e-05
Newton iteration 0: r (abs) = 2.199e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.946e+19 (tol = 1.000e-10) r (rel) = 1.340e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.765e+16 (tol = 1.000e-10) r (rel) = 1.257e+00 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 2.421e+10 (tol = 1.000e-10) r (rel) = 1.101e-06 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.642e+08 (tol = 1.000e-10) r (rel) = 7.470e-09 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Step:  15 // t = 3.0000000000000004e-05
Newton iteration 0: r (abs) = 2.042e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.500e+19 (tol = 1.000e-10) r (rel) = 1.224e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.037e+16 (tol = 1.000e-10) r (rel) = 9.977e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.350e+10 (tol = 1.000e-10) r (rel) = 6.611e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  16 // t = 3.2000000000000005e-05
Newton iteration 0: r (abs) = 1.887e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.133e+19 (tol = 1.000e-10) r (rel) = 1.130e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.515e+16 (tol = 1.000e-10) r (rel) = 8.028e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 7.651e+09 (tol = 1.000e-10) r (rel) = 4.055e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  17 // t = 3.4000000000000007e-05
Newton iteration 0: r (abs) = 1.734e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.825e+19 (tol = 1.000e-10) r (rel) = 1.053e+03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.131e+16 (tol = 1.000e-10) r (rel) = 6.522e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 4.363e+09 (tol = 1.000e-10) r (rel) = 2.516e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  18 // t = 3.600000000000001e-05
Newton iteration 0: r (abs) = 1.582e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.561e+19 (tol = 1.000e-10) r (rel) = 9.866e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 8.425e+15 (tol = 1.000e-10) r (rel) = 5.324e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 2.469e+09 (tol = 1.000e-10) r (rel) = 1.560e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  19 // t = 3.800000000000001e-05
Newton iteration 0: r (abs) = 1.432e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.331e+19 (tol = 1.000e-10) r (rel) = 9.291e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 6.221e+15 (tol = 1.000e-10) r (rel) = 4.343e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.385e+09 (tol = 1.000e-10) r (rel) = 9.671e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  20 // t = 4.000000000000001e-05
Newton iteration 0: r (abs) = 1.284e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.126e+19 (tol = 1.000e-10) r (rel) = 8.771e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 4.520e+15 (tol = 1.000e-10) r (rel) = 3.521e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 7.515e+08 (tol = 1.000e-10) r (rel) = 5.854e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  21 // t = 4.2000000000000004e-05
Newton iteration 0: r (abs) = 1.138e+16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 9.423e+18 (tol = 1.000e-10) r (rel) = 8.279e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 3.205e+15 (tol = 1.000e-10) r (rel) = 2.816e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 4.110e+08 (tol = 1.000e-10) r (rel) = 3.611e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  22 // t = 4.4000000000000006e-05
Newton iteration 0: r (abs) = 9.962e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 7.762e+18 (tol = 1.000e-10) r (rel) = 7.792e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.198e+15 (tol = 1.000e-10) r (rel) = 2.207e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 2.438e+08 (tol = 1.000e-10) r (rel) = 2.447e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  23 // t = 4.600000000000001e-05
Newton iteration 0: r (abs) = 8.594e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 6.262e+18 (tol = 1.000e-10) r (rel) = 7.287e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.444e+15 (tol = 1.000e-10) r (rel) = 1.680e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.865e+08 (tol = 1.000e-10) r (rel) = 2.170e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  24 // t = 4.800000000000001e-05
Newton iteration 0: r (abs) = 7.295e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 4.923e+18 (tol = 1.000e-10) r (rel) = 6.748e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 8.980e+14 (tol = 1.000e-10) r (rel) = 1.231e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.718e+08 (tol = 1.000e-10) r (rel) = 2.355e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  25 // t = 5.000000000000001e-05
Newton iteration 0: r (abs) = 6.080e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 3.750e+18 (tol = 1.000e-10) r (rel) = 6.167e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 5.232e+14 (tol = 1.000e-10) r (rel) = 8.605e-02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.685e+08 (tol = 1.000e-10) r (rel) = 2.771e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  26 // t = 5.200000000000001e-05
Newton iteration 0: r (abs) = 4.968e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.752e+18 (tol = 1.000e-10) r (rel) = 5.539e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.824e+14 (tol = 1.000e-10) r (rel) = 5.684e-02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.661e+08 (tol = 1.000e-10) r (rel) = 3.344e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  27 // t = 5.400000000000001e-05
Newton iteration 0: r (abs) = 3.971e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.935e+18 (tol = 1.000e-10) r (rel) = 4.873e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.395e+14 (tol = 1.000e-10) r (rel) = 3.514e-02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.699e+08 (tol = 1.000e-10) r (rel) = 4.279e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  28 // t = 5.6000000000000006e-05
Newton iteration 0: r (abs) = 3.100e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.297e+18 (tol = 1.000e-10) r (rel) = 4.184e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 6.249e+13 (tol = 1.000e-10) r (rel) = 2.016e-02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.725e+08 (tol = 1.000e-10) r (rel) = 5.567e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  29 // t = 5.800000000000001e-05
Newton iteration 0: r (abs) = 2.359e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 8.245e+17 (tol = 1.000e-10) r (rel) = 3.494e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.512e+13 (tol = 1.000e-10) r (rel) = 1.065e-02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.662e+08 (tol = 1.000e-10) r (rel) = 7.042e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  30 // t = 6.000000000000001e-05
Newton iteration 0: r (abs) = 1.749e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 4.952e+17 (tol = 1.000e-10) r (rel) = 2.832e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 8.988e+12 (tol = 1.000e-10) r (rel) = 5.140e-03 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.521e+08 (tol = 1.000e-10) r (rel) = 8.700e-08 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  31 // t = 6.200000000000002e-05
Newton iteration 0: r (abs) = 1.261e+15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.801e+17 (tol = 1.000e-10) r (rel) = 2.222e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.846e+12 (tol = 1.000e-10) r (rel) = 2.257e-03 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.591e+08 (tol = 1.000e-10) r (rel) = 1.262e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  32 // t = 6.400000000000001e-05
Newton iteration 0: r (abs) = 8.835e+14 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.488e+17 (tol = 1.000e-10) r (rel) = 1.685e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 7.938e+11 (tol = 1.000e-10) r (rel) = 8.984e-04 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.567e+08 (tol = 1.000e-10) r (rel) = 1.774e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  33 // t = 6.6e-05
Newton iteration 0: r (abs) = 6.015e+14 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 7.421e+16 (tol = 1.000e-10) r (rel) = 1.234e+02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.945e+11 (tol = 1.000e-10) r (rel) = 3.234e-04 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.625e+08 (tol = 1.000e-10) r (rel) = 2.702e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  34 // t = 6.800000000000001e-05
Newton iteration 0: r (abs) = 3.977e+14 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 3.469e+16 (tol = 1.000e-10) r (rel) = 8.723e+01 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 4.184e+10 (tol = 1.000e-10) r (rel) = 1.052e-04 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.595e+08 (tol = 1.000e-10) r (rel) = 4.011e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  35 // t = 7.000000000000001e-05
Newton iteration 0: r (abs) = 2.554e+14 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.521e+16 (tol = 1.000e-10) r (rel) = 5.954e+01 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 7.911e+09 (tol = 1.000e-10) r (rel) = 3.097e-05 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.544e+08 (tol = 1.000e-10) r (rel) = 6.044e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  36 // t = 7.200000000000002e-05
Newton iteration 0: r (abs) = 1.594e+14 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 6.257e+15 (tol = 1.000e-10) r (rel) = 3.925e+01 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.323e+09 (tol = 1.000e-10) r (rel) = 8.298e-06 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.591e+08 (tol = 1.000e-10) r (rel) = 9.981e-07 (tol = 1.000e-06)
Newton solver finished in 3 iterations and 3 linear solver iterations.
Step:  37 // t = 7.400000000000001e-05
Newton iteration 0: r (abs) = 9.667e+13 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.419e+15 (tol = 1.000e-10) r (rel) = 2.502e+01 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 2.558e+08 (tol = 1.000e-10) r (rel) = 2.646e-06 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.565e+08 (tol = 1.000e-10) r (rel) = 1.619e-06 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.583e+08 (tol = 1.000e-10) r (rel) = 1.637e-06 (tol = 1.000e-06)
Newton iteration 5: r (abs) = 1.572e+08 (tol = 1.000e-10) r (rel) = 1.626e-06 (tol = 1.000e-06)
Newton iteration 6: r (abs) = 1.541e+08 (tol = 1.000e-10) r (rel) = 1.594e-06 (tol = 1.000e-06)
Newton iteration 7: r (abs) = 1.530e+08 (tol = 1.000e-10) r (rel) = 1.583e-06 (tol = 1.000e-06)
Newton iteration 8: r (abs) = 1.599e+08 (tol = 1.000e-10) r (rel) = 1.654e-06 (tol = 1.000e-06)
Newton iteration 9: r (abs) = 1.562e+08 (tol = 1.000e-10) r (rel) = 1.616e-06 (tol = 1.000e-06)
Newton iteration 10: r (abs) = 1.535e+08 (tol = 1.000e-10) r (rel) = 1.588e-06 (tol = 1.000e-06)
Newton iteration 11: r (abs) = 1.556e+08 (tol = 1.000e-10) r (rel) = 1.610e-06 (tol = 1.000e-06)
Newton iteration 12: r (abs) = 1.585e+08 (tol = 1.000e-10) r (rel) = 1.639e-06 (tol = 1.000e-06)
Newton iteration 13: r (abs) = 1.647e+08 (tol = 1.000e-10) r (rel) = 1.704e-06 (tol = 1.000e-06)
Newton iteration 14: r (abs) = 1.599e+08 (tol = 1.000e-10) r (rel) = 1.654e-06 (tol = 1.000e-06)
Newton iteration 15: r (abs) = 1.558e+08 (tol = 1.000e-10) r (rel) = 1.612e-06 (tol = 1.000e-06)
Newton iteration 16: r (abs) = 1.498e+08 (tol = 1.000e-10) r (rel) = 1.549e-06 (tol = 1.000e-06)
Newton iteration 17: r (abs) = 1.537e+08 (tol = 1.000e-10) r (rel) = 1.589e-06 (tol = 1.000e-06)
Newton iteration 18: r (abs) = 1.558e+08 (tol = 1.000e-10) r (rel) = 1.612e-06 (tol = 1.000e-06)
Newton iteration 19: r (abs) = 1.559e+08 (tol = 1.000e-10) r (rel) = 1.613e-06 (tol = 1.000e-06)
Newton iteration 20: r (abs) = 1.562e+08 (tol = 1.000e-10) r (rel) = 1.616e-06 (tol = 1.000e-06)
Newton iteration 21: r (abs) = 1.526e+08 (tol = 1.000e-10) r (rel) = 1.579e-06 (tol = 1.000e-06)
Newton iteration 22: r (abs) = 1.562e+08 (tol = 1.000e-10) r (rel) = 1.616e-06 (tol = 1.000e-06)
Newton iteration 23: r (abs) = 1.618e+08 (tol = 1.000e-10) r (rel) = 1.673e-06 (tol = 1.000e-06)
Newton iteration 24: r (abs) = 1.602e+08 (tol = 1.000e-10) r (rel) = 1.657e-06 (tol = 1.000e-06)
Newton iteration 25: r (abs) = 1.580e+08 (tol = 1.000e-10) r (rel) = 1.634e-06 (tol = 1.000e-06)
Newton iteration 26: r (abs) = 1.518e+08 (tol = 1.000e-10) r (rel) = 1.570e-06 (tol = 1.000e-06)
Newton iteration 27: r (abs) = 1.508e+08 (tol = 1.000e-10) r (rel) = 1.560e-06 (tol = 1.000e-06)
Newton iteration 28: r (abs) = 1.541e+08 (tol = 1.000e-10) r (rel) = 1.594e-06 (tol = 1.000e-06)
Newton iteration 29: r (abs) = 1.625e+08 (tol = 1.000e-10) r (rel) = 1.681e-06 (tol = 1.000e-06)
Newton iteration 30: r (abs) = 1.618e+08 (tol = 1.000e-10) r (rel) = 1.673e-06 (tol = 1.000e-06)
Newton iteration 31: r (abs) = 1.591e+08 (tol = 1.000e-10) r (rel) = 1.645e-06 (tol = 1.000e-06)
Newton iteration 32: r (abs) = 1.551e+08 (tol = 1.000e-10) r (rel) = 1.605e-06 (tol = 1.000e-06)
Newton iteration 33: r (abs) = 1.603e+08 (tol = 1.000e-10) r (rel) = 1.658e-06 (tol = 1.000e-06)
Newton iteration 34: r (abs) = 1.561e+08 (tol = 1.000e-10) r (rel) = 1.614e-06 (tol = 1.000e-06)
Newton iteration 35: r (abs) = 1.551e+08 (tol = 1.000e-10) r (rel) = 1.605e-06 (tol = 1.000e-06)
Newton iteration 36: r (abs) = 1.591e+08 (tol = 1.000e-10) r (rel) = 1.646e-06 (tol = 1.000e-06)
Newton iteration 37: r (abs) = 1.555e+08 (tol = 1.000e-10) r (rel) = 1.609e-06 (tol = 1.000e-06)
Newton iteration 38: r (abs) = 1.538e+08 (tol = 1.000e-10) r (rel) = 1.591e-06 (tol = 1.000e-06)
Newton iteration 39: r (abs) = 1.588e+08 (tol = 1.000e-10) r (rel) = 1.643e-06 (tol = 1.000e-06)
Newton iteration 40: r (abs) = 1.519e+08 (tol = 1.000e-10) r (rel) = 1.572e-06 (tol = 1.000e-06)
Newton iteration 41: r (abs) = 1.484e+08 (tol = 1.000e-10) r (rel) = 1.535e-06 (tol = 1.000e-06)
Newton iteration 42: r (abs) = 1.538e+08 (tol = 1.000e-10) r (rel) = 1.591e-06 (tol = 1.000e-06)
Newton iteration 43: r (abs) = 1.594e+08 (tol = 1.000e-10) r (rel) = 1.649e-06 (tol = 1.000e-06)
Newton iteration 44: r (abs) = 1.590e+08 (tol = 1.000e-10) r (rel) = 1.645e-06 (tol = 1.000e-06)
Newton iteration 45: r (abs) = 1.548e+08 (tol = 1.000e-10) r (rel) = 1.601e-06 (tol = 1.000e-06)
Newton iteration 46: r (abs) = 1.525e+08 (tol = 1.000e-10) r (rel) = 1.578e-06 (tol = 1.000e-06)
Newton iteration 47: r (abs) = 1.560e+08 (tol = 1.000e-10) r (rel) = 1.614e-06 (tol = 1.000e-06)
Newton iteration 48: r (abs) = 1.574e+08 (tol = 1.000e-10) r (rel) = 1.628e-06 (tol = 1.000e-06)
Newton iteration 49: r (abs) = 1.531e+08 (tol = 1.000e-10) r (rel) = 1.583e-06 (tol = 1.000e-06)
Newton iteration 50: r (abs) = 1.564e+08 (tol = 1.000e-10) r (rel) = 1.618e-06 (tol = 1.000e-06)
Traceback (most recent call last):
File "mwe.py", line 134, in <module>
solver.solve(Problem(), u_vec.vector())
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
***
***
*** 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
*** -------------------------------------------------------------------------

1 Like

The solution looks pretty much as I expect until some point in time when the solver does not converge.

Strange, I just copied and pasted your code, ran with 4 processes, and it worked…

• In the first timestep the solver takes much longer to converge than in the other ones (18 iterations vs 3 or 4 in the subsequent ones).

This is because your initial state/guess is probably very poor. Subsequent initial guesses are using the solution from the previous time step, which should be deeper in an attractor region.

• In all timeteps in the first iteration of the solver the residual increases before decreasing and converging to the solution I don’t know if this is common or if I should check if my residual is well defined.

Is your chosen value of dt_n appropriate?

Newton iteration 0: r (abs) = 3.977e+14 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)

Your residuals’ absolute norms are huge. I’d highly recommend you non dimensionalise for n0 and u0.

There are some general tips and tricks for solving nonlinear problems I wrote a while ago here.

I don’t really know what else to suggest. But if it were me, I’d manufacture a solution such that I could check convergence rates. That way I’d have far more confidence that I’m solving the equations I think I’m solving.

1 Like

True, I am sorry, this only happens if you increase the final time.

Alright, I will take that into account and see if I can solve it.

I am not sure to be honest I will try different time steps and see how that behaves.

YesI thought about it but I didn’t have time to implement it yet, it’s definitely necessary. I’ll see how the problem behaves after non dimensionalisation.

Great thanks for the tips, I’ll implement manufactured solutions to check everything is working fine.

1 Like