Newton solver not converging for simple heat problem

General
1

Hi – I have a time-dependent heat equation that I am trying to solve, but Newton solver was not converging, so I tried solving a simpler problem (f=1 and Gaussian initial condition), but it is still not converging, which makes me think I have set something up wrong. The domain is just a cube (which I think I have set up right?), and I have homogeneous Dirichlet b.c. Even if I am extremely liberal with the relative tolerance, it still doesn’t converge. Below is a minimal working example. Any help is greatly appreciated!

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh, io, plot, nls, log
log.set_log_level(log.LogLevel.INFO)
import ufl

# Define initial condition
def initial_condition(x):
    return np.exp(-(x[0]**2+x[1]**2+x[2]**2))

# Parameter for reaction-diffusion
D = 800 # µm^2/s

def rxn_dif(d,rho):
    
    # Define temporal parameters
    t = 0 # Start time
    T = 10.0 # Final time
    num_steps = 500     
    dt = T / num_steps # time step size

    # MESH

    nx, ny, nz = 50, 50, 50
    domain = mesh.create_box(MPI.COMM_WORLD, [np.array([-10,-10,-10]), np.array([10, 10, 10])], 
                                   [nx, ny, nz], mesh.CellType.tetrahedron)
    
    V = fem.FunctionSpace(domain, ("CG", 2))


    # ------------------------------------------------------------------------------------------------

    # INITIAL CONDITION

    u_n = fem.Function(V)
    u_n.name = "u_n"
    u_n.interpolate(initial_condition)

    # ------------------------------------------------------------------------------------------------

    # BOUNDARY CONDITIONS

    fdim = domain.topology.dim - 1
    boundary_facets = mesh.locate_entities_boundary(
        domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
    bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)

    # ------------------------------------------------------------------------------------------------

    # Define solution variable and interpolate initial solution

    uh = fem.Function(V)
    uh.name = "uh"
    uh.interpolate(initial_condition)

    # ------------------------------------------------------------------------------------------------

    # VARIATIONAL FORM
    
    x = ufl.SpatialCoordinate(domain)
    v = ufl.TestFunction(V)
    f = 1.0
    F = uh * v * ufl.dx + dt*ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - (u_n + dt * f) * v * ufl.dx
    problem = fem.petsc.NonlinearProblem(F, uh, bcs=[bc])

    # ------------------------------------------------------------------------------------------------

    # SOLVE

    solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
    
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-4
    solver.report = True

    for i in range(num_steps):
        t += dt
        print(t)

        # Solve linear problem
        n, converged = solver.solve(uh)
        assert(converged)
        print(f"Number of interations: {n:d}")
        uh.x.scatter_forward()

        # Update solution at previous time step (u_n)
        u_n.x.array[:] = uh.x.array
        
    return uh
    

print(rxn_dif(5,1))
2

Your code seems to work fine for me. Could you add

    ...
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-4
    solver.report = True

    # -- Add this
    pfx = solver.krylov_solver.getOptionsPrefix()
    opt = PETSc.Options()
    opt[f"{pfx}ksp_view"] = None
    solver.krylov_solver.setFromOptions()
    # ---

    for i in range(num_steps):
    ...

And let us know what solver is being used.

3

Hi, thanks for the response. It is using the Krylov solver, if I add the code that you indicated. Output is:

2023-01-23 09:58:03.767 ( 4.955s) [main ] petsc.cpp:677 INFO| PETSc Krylov solver starting to solve system.

It does not advance beyond the first timestep.

4

Very strange. So you don’t even see the output from the PETSc KSP view? E.g. I see the default ksp_type as preonly and pc_type as lu.

5

No, unfortunately I only see the print(t) output (0.02) and then the output I mentioned in the last comment. I have gotten the heat equation code from the tutorial to work, so I think it is some issue with the nonlinear solver on my system.