I am working on solving the following PDE: $$\left(\mu_{x}\frac{\partial^{2}u}{\partial x^{2}}+\mu_{y}\frac{\partial^{2}u}{\partial y^{2}}\right)=f(x,y) \tag 1$$

Which is then discretised:

$$- \mu_{x} \left(\frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{h^2}\right) - \mu_{y} \left(\frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{h^2}\right) = f_{i,j}$$

once all the nonlinear functions have been computed at the current guess $u$, we have $A$ as the Jacobian matrix in newtons method and $F$ for the right-hand side.

We now have a linear equation system $JU = F$. I am now trying to compute $U$ using sp.sparse.linalg.cg however, when I run it, it seems to give me 0 iterations and a result of 0. 0 usually indicates convergence but I can't see how it's managed to do this in 0 iterations. Would appreciate any insight into where I might be doing wrong with this.

I've provided a full reproducible example: solver.py

import scipy as sp
import numpy as np
from scipy.sparse.linalg import norm

def NewtonSys(fnon, jac, x0, tol, maxk, *fnonargs):
    # fnon     - name of the nonlinear function f(x)
    # jac      - name of the Jacobian function J(x)
    # x0       - initial guess for the solution x0
    # tol      - stopping tolerance for Newton's iteration
    # maxk     - maximum number of Newton iterations before stopping
    # fnonargs - optional arguments that will be passed to the nonlinear function (useful for additional function parameters)

    k = 0
    x = x0

    F = eval(fnon)(x,*fnonargs)

    Fsize = F.get_shape()
    n = Fsize[0]

    normF = (F.data ** 2).sum()

    print(' k    f(xk)')

    # Main Newton loop
    while (normF > tol and k <= maxk):
        # Evaluate Jacobian matrix
        J = eval(jac)(x,n,fnon,F,*fnonargs)

        # Take Newton step by solving the tangent problem
        delta = np.linalg.solve(J.toarray(),F.toarray())
        x = x - delta

        F = eval(fnon)(x,*fnonargs) 
        normF = (F.data ** 2).sum()

        print('{0:2.0f}  {1:2.2e}'.format(k, normF))
        k += 1

    if (k >= maxk):
        print('Not converged')
        return x, J, F

def matrix( mu_x, mu_y, n ):

  N = n**2

  # We use a list-of-lists format for more efficient assembly
  A = sp.sparse.lil_array((N, N), dtype=np.float64)
  h = 1. / (n-1);

  stencil = np.array([2 * mu_x + 2 * mu_y, mu_x * -1., mu_x * -1., mu_y * -1., mu_y * -1.]) / h**2

  # Loop over each internal node in the grid, i,j = 1,2,...,n-2
  # starts at the second row to avoid boundary and ends at n-2 to avoid boundary as well
  for i in range(1,n-1):
    for j in range(1,n-1):
      # Find k-indices of the four neighbouring nodes
      localStencilIndices = np.array([indexFD(i,j,n), indexFD(i+1,j,n), indexFD(i-1,j,n), indexFD(i,j+1,n), indexFD(i,j-1,n)])

      # Add the local stencil for node (x_i,y_j) to the matrix
      currentRow = indexFD(i,j,n)
      for m in range (0,5):
        A[currentRow, localStencilIndices[m]] = A[currentRow, localStencilIndices[m]] + stencil[m]

  # After the matrix A has been assembled, we convert it to the column-major format
  # for more efficient computations
  A = A.tocsc()

  return (A)

def indexFD(i, j, n):
  return ( i*n + j )

import numpy as np

def boundaryConditions(n):

  extNodes = []

  # Find nodes at the boundary of the square
  for j in range(0,n):
    extNodes.append( indexFD(j,0,n) )
    extNodes.append( indexFD(j,n-1,n) )
    extNodes.append( indexFD(0,j,n) )
    extNodes.append( indexFD(n-1,j,n) )

  extNodes = np.unique(extNodes)
  intNodes = np.setdiff1d(np.arange(0,n**2), extNodes);

  return intNodes, extNodes

import scipy as sp
import copy

def fdJacobian(x,n,fnon,F0,*fnonargs):
  # We replace the dense NumPy array with a sparse SciPy matrix
  # (list-of-lists format is used for assembly)

  J = sp.sparse.lil_array((n, n), dtype=np.float64)
  h = 10e-8

  for k in range(0,n):
    xb = copy.deepcopy(x)
    xb[k,0] = xb[k,0] + h

    F = eval(fnon)(xb,*fnonargs)

    for i in range(0,n):
      J[i,k] = (F[i,0] - F0[i,0]) / h

  # Return the Jacobian in column-major sparse format
  return J.tocsc()

import math
def source_function(x, y, h):
      x = x * h
      y = y * h
      if x >= 0.1 and x <= 0.3 and y >= 0.1 and y <= 0.3:
        return 1
        return 0

def solve_rhs( u, A, n ):

  N = n**2
  capitalF = sp.sparse.lil_array((N, 1), dtype=np.float64)

  # Find the internal nodes and the boundary nodes
  intNodes, extNodes  = boundaryConditions(n)

  Au = A[intNodes,:] @ u # Sparse matrix-vector multiplication
  littleF = sp.sparse.lil_array((N, 1), dtype=np.float64)

  for k in intNodes:
      i, j = index_to_coordinates(k,n)
      f_ij = source_function(i, j, 0.1)
      littleF[k] = f_ij

  capitalF[intNodes] = Au - littleF[intNodes]

  # Set value of U at boundary nodes (this imposes u=0 at the boundary)
  capitalF[extNodes] = u[extNodes]

  return capitalF.tocsc()

def Solver(mu_x, mu_y, h):

  n = int(1./h + 1)  # dimension of spatial mesh in each dimension
  N = n**2      # total number of grid points

  # Initial guess for Newton is the zero vector
  u0 = sp.sparse.csc_array( (N,1), dtype=np.float64 )

  # Assemble the Laplacian matrix
  A = matrix(mu_x, mu_y, n)

  # We plot the sparsity pattern of the matrix A to check that it has been
  # correctly assembled.
  fig = plt.figure()
  ax = fig.add_subplot(111)
  cax = ax.matshow(A.toarray(), vmin=-4./h*mu_x, vmax=4./h*mu_y, cmap='coolwarm')
  plt.title("Sparsity pattern")

  # Newton iteration to solve nonlinear PDE
  u, J, F = NewtonSys("solve_rhs", "fdJacobian", u0, 1e-6, 100, A, n )

  return u, J, F

h = 0.1
#mu_x = 1.0, mu_y = 1.0
sol, J, F = Solver(1.0, 1.0, h)

# Solve AU = F 
tolerance = 1e-6
max_iterations = 1000

def callback(xk):
    callback.num_iterations += 1

callback.num_iterations = 0
x1, info1 = sp.sparse.linalg.cg(J, F.todense(), tol=tolerance, maxiter=max_iterations, callback=callback)
print("Iterations: ", callback.num_iterations)
print("Info: ", info1)
  • $\begingroup$ why there is a newton step solving nonlinear PDE (see your comment)? Isn't just getting the matrix A and solve by cg? $\endgroup$
    – Yimin
    Dec 11, 2023 at 22:58
  • $\begingroup$ There is a lot going on here, but it looks like you have both a finite difference function and a sparse array to compute the Laplacian being passed to a Newton solver for a linear problem. Try to reduce this to a minimal working example and I suspect you will figure out your problems. $\endgroup$
    – whpowell96
    Dec 12, 2023 at 17:04


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.