Cannot run Fenics code with fine mesh (memory issue?)

[mekong]

Dear all,
I am using a Fenics code in Python, with Docker, which solves a PDE with mixed finite-elements. The mesh is given by the region between two disks, both centered at the origin, with radii r=1 and R=2,. The code runs with no problem if the mesh resolution is >= 0.05, with the following output

$python3 example.py
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.002e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 1.805e-01 (tol = 1.000e-10) r (rel) = 1.802e-03 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 8.813e-05 (tol = 1.000e-10) r (rel) = 8.799e-07 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 4.045e-11 (tol = 1.000e-10) r (rel) = 4.038e-13 (tol = 1.000e-09)
  Newton solver finished in 3 iterations and 3 linear solver iterations.

However, if I set the resolution to r=0.025 or smaller, I get the following error

$python3 example.py
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.917e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

Traceback (most recent call last):
  File "example.py", line 206, 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 successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Do you have any idea of what is going on ?
I suspect that this is a memory issue, because if I run $top while the code is running, I get that the code is using 8 GB of RAM

PID    COMMAND      %CPU  TIME     #TH   #WQ  #PORT MEM    PURG  CMPRS  PGRP  PPID  STATE    BOOSTS                %CPU_ME %CPU_OTHRS UID  FAULTS    COW      MSGSENT     MSGRECV     SYSBSD      SYSMACH
95703  com.apple.Vi 100.4 03:15:43 19/1  1    78    8065M  0B    9687M  95703 1     running  *5[3]                 0.00000 0.00000    501  9318288   125      934         206         37862101+   751

If so, how may I increase the memory of a docker container ? I tried docker update --kernel-memory 16G my_fenics_container but nothing changes.

Thank you

[dokken]

@mekong

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

The error is quite clear.
Please search the forum for this message.

UMFPACK error: out of memory despite system having free memory

Hey, this could indeed be a bug of umfpack. As far as I know this is because fenics is using a 32 bit version of this (someone please correct me if I’m wrong). What you can do, however, is to use mumps, which is 64 bit, and can support more memory. For a bilinear form a, linear form L and boundary conditions given in bcs, the solve command for this looks like solve(a==L, u, bcs, solver_parameters={'linear_solver' : 'mumps'}) If you want more freedom in configuring mumps, I’ll refer you to my…

Large Mesh Error PETSc Error 76 - Adjoint Calculation

So this error message is related to the UMFPACK bug i referenced above. However, I notice that you get an actual error when running: solve(A, w.vector(), L, "mumps") Temporary fix would be to change this to: solver_lin = LUSolver("mumps") solver_lin.solve(A, w.vector(), L) Edit: PR added

[mekong]

@dokken

The error is quite clear.

Oh sorry I had not seen that! I will look into the links that you provided, thank you!

[mekong]

I tried to use the suggestion given in the first link by using mumps as a liner solver

params = {'nonlinear_solver': 'newton',
          'newton_solver':
              {
                  'linear_solver': 'mumps',
              }
          }
solver.parameters.update( params )
solver.solve()

but it does not work.
I also tried this suggestion but I could not even install pyadjoint

$  sudo pip3 install git+https://github.com/dolfin-adjoint/pyadjoint.git
WARNING: The directory '/home/fenics/.cache/pip/http' or its parent directory is not owned by the current user and the cache has been disabled. Please check the permissions and owner of that directory. If executing pip with sudo, you may want sudo's -H flag.
WARNING: The directory '/home/fenics/.cache/pip' or its parent directory is not owned by the current user and caching wheels has been disabled. check the permissions and owner of that directory. If executing pip with sudo, you may want sudo's -H flag.
Collecting git+https://github.com/dolfin-adjoint/pyadjoint.git
  Cloning https://github.com/dolfin-adjoint/pyadjoint.git to /tmp/pip-req-build-_1bxm1ih
  Running command git clone -q https://github.com/dolfin-adjoint/pyadjoint.git /tmp/pip-req-build-_1bxm1ih
Collecting scipy>=1.0
  Downloading https://files.pythonhosted.org/packages/c8/89/63171228d5ced148f5ced50305c89e8576ffc695a90b58fe5bb602b910c2/scipy-1.5.4-cp36-cp36m-manylinux1_x86_64.whl (25.9MB)
     |████████████████████████████████| 25.9MB 16.0MB/s 
ERROR: Could not find a version that satisfies the requirement checkpoint_schedules (from pyadjoint-ad==2023.0.0) (from versions: none)
ERROR: No matching distribution found for checkpoint_schedules (from pyadjoint-ad==2023.0.0)

Also, using the latest Fenics version does not help either: docker run -ti -v $(pwd):/home/fenics/shared -w /home/fenics/shared/ --rm quay.io/fenicsproject/dev:latest gives, when running,

  Process 0: Newton iteration 0: r (abs) = 3.626e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
  File "example.py", line 206, 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 linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 3
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------

Any suggestion ?

[dokken]

First of all. When you say it doesn’t work, what is the traceback?

Secondly, no need to install pyadjoint, as it has nothing to do with the problem you are trying to solve.

[mekong]

@dokken

First of all. When you say it doesn’t work, what is the traceback?

Here it is

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.917e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
  File "example.py", line 230, 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 linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

[dokken]

Well, now you are getting another error, not out of memory.

@mekong

*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).

You can again search for this error on the forum ,getting results like:

ArityMismatch error

Thanks for the insightful take. I appreciate it. Kindly provide your mesh file (in this case .xml file) with a minimum working example so that I can reproduce the error and could possibly give any suggestion. According to @dokken in this post PETSc Krylov Solver RuntimeError with fine mesh , this information means that the residual is zero at the first iteration, which usually indicates that the initial guess solves your problem. This was also reported for another case at the old Q&A forum. …

PETSc Krylov Solver RuntimeError with fine mesh

This means that the residual is zero at the first iteration, which usually indicates that the initial guess solves your problem. This was also reported for another case at the old Q&A forum.

[mekong]

Following the suggestions in the first link,
I used

params = {'nonlinear_solver': 'newton', 'newton_solver': { 'linear_solver' : 'umfpack',  },}

and we are back to an out-of-memory error:

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.917e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

Traceback (most recent call last):
  File "example.py", line 230, 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 successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

I tried replacing umfpack with lu and I get the error

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.917e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

Traceback (most recent call last):
  File "example.py", line 230, 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 successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

[dokken]

@mekong

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.917e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

As you can see “lu”-> “umfpack”, so you need to use “mumps” or “superlu”.

As you haven’t provided any code, I cannot give you any further help.

[mekong]

@dokken

so you need to use “mumps” or “superlu”.

I tried mumps above and it gave the error I reported . With superlu I get

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.917e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Killed

[mekong]

@dokken

As you haven’t provided any code, I cannot give you any further help.

I made a minimal working example for you, mwe.py:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import argparse
import ufl as ufl

#### function definitions
i, j, k, l = ufl.indices(4)

def my_norm(x):
    return (sqrt(np.dot(x, x)))

def ufl_norm(x):
    return(sqrt(ufl.dot(x, x)))

def X(z):
    x = ufl.SpatialCoordinate(mesh)
    return as_tensor([x[0], x[1], z])

def e(omega):
    return as_tensor([[1, 0, omega[0]], [0, 1, omega[1]]])

def normal(omega):
    return as_tensor(cross(e(omega)[0], e(omega)[1]) /  ufl_norm(cross(e(omega)[0], e(omega)[1])) )

def b(omega):
    return as_tensor((normal(omega))[k] * (e(omega)[i, k]).dx(j), (i,j))

def g(omega):
    return as_tensor([[1+ (omega[0])**2, (omega[0])*(omega[1])],[(omega[0])*(omega[1]), 1+ (omega[1])**2]])

def g_c(omega):
    return ufl.inv(g(omega))

def detg(omega):
    return ufl.det(g(omega))

def abs_detg(omega):
    return np.abs(ufl.det(g(omega)))

def sqrt_detg(omega):
    return sqrt(detg(omega))

def sqrt_abs_detg(omega):
    return sqrt(abs_detg(omega))

def w():
    x = ufl.SpatialCoordinate(mesh)
    return as_tensor([-x[1], x[0]])

def sqrt_deth(omega):
    return(sqrt((w())[i]*(w())[j]*g(omega)[i, j]))

def calc_normal_cg2(mesh):
    n = FacetNormal(mesh)
    V = VectorFunctionSpace(mesh, "CG", 2)
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(u, v) * ds
    l = inner(n, v) * ds
    A = assemble(a, keep_diagonal=True)
    L = assemble(l)

    A.ident_zeros()
    nh = Function(V)
    solve(A, nh.vector(), L)
    return nh

def n(omega):
    u = n_overline
    return as_tensor(u[k]/sqrt(g(omega)[i,j]*u[i]*u[j]), (k))

def H(omega):
    return (0.5 * g_c(omega)[i, j]*b(omega)[j, i])

def K(omega):
    return(ufl.det(as_tensor(b(omega)[i,k]*g_c(omega)[k,j], (i, j))))

def Gamma(omega):
    return as_tensor(0.5 * g_c(omega)[i,l] * ( (g(omega)[l, k]).dx(j) + (g(omega)[j, l]).dx(k) - (g(omega)[j, k]).dx(l) ), (i, j, k))

def Nabla_v(u, omega):
    return as_tensor((u[i]).dx(j) + u[k] * Gamma(omega)[i, k, j], (i, j))

def Nabla_f(f, omega):
    return as_tensor((f[i]).dx(j) - f[k] * Gamma(omega)[k, i, j], (i, j))
####


parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
args = parser.parse_args()

#read mesh
mesh=Mesh()
with XDMFFile((args.input_directory) + "/triangle_mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile((args.input_directory) + "/line_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")

n_overline = FacetNormal( mesh )

# Define boundaries and obstacle
boundary = 'on_boundary'
boundary_r = 'on_boundary && sqrt(pow(x[0], 2) + pow(x[1], 2)) < (1.0+2.0)/2.0'
boundary_R = 'on_boundary && sqrt(pow(x[0], 2) + pow(x[1], 2)) > (1.0+2.0)/2.0'


r = 1.0
R = 2.0
kappa = 1.0
alpha = 1000.0
sigma0 = 1.0
C = 0.5

set_log_level(20)


# Define function spaces
P_z = FiniteElement('P', triangle, 1)
P_omega = VectorElement('P', triangle, 3)
element = MixedElement([P_z, P_omega])
Q_z_omega = FunctionSpace(mesh, element)
Q_z = Q_z_omega.sub(0).collapse()
Q_omega = Q_z_omega.sub(1).collapse()

# Create XDMF files for visualization output
xdmffile_z = XDMFFile('z.xdmf')
xdmffile_omega = XDMFFile('omega.xdmf')

#read an object with label subdomain_id from xdmf file and assign to it the ds `ds_inner`
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_r = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_R = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)

# Define trial and test functions
nu_z, nu_omega = TestFunctions(Q_z_omega)
z_omega = Function(Q_z_omega)
sigma = Function(Q_z)
z, omega = split(z_omega)


class grad_r_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = C * x[0] / my_norm( x )
        values[1] = C * x[1] / my_norm( x )
        return (2,)

class grad_R_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = 0
        values[1] = 0
    def value_shape(self):
        return (2,)

class sigma_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = sigma0
    def value_shape(self):
        return (1,)


# Define expressions used in variational forms
kappa = Constant(kappa)
sigma.interpolate(sigma_Expression(element=Q_z.ufl_element()))
grad_r = interpolate(grad_r_Expression(element=Q_omega.ufl_element()), Q_omega)
grad_R = interpolate(grad_R_Expression(element=Q_omega.ufl_element()), Q_omega)

# Define functional for variational problem
F_z = ( kappa * ( g_c(omega)[i, j] * (H(omega).dx(j)) * (nu_z.dx(i)) - 2.0 * H(omega) * ( (H(omega))**2 - K(omega) ) * nu_z ) + sigma * H(omega) * nu_z ) * sqrt_detg(omega) * dx \
    - ( kappa * (n(omega))[i] * nu_z * (H(omega).dx(i)) ) * sqrt_deth(omega) * ds
F_omega = ( - z * Nabla_v(nu_omega, omega)[i, i] - omega[i] * nu_omega[i] ) *  sqrt_detg(omega) * dx + \
          ( (n(omega))[i] * g(omega)[i, j] * z * nu_omega[j] ) * sqrt_deth(omega) * ds
#Nitche's part of the functional
F_N = alpha * ( \
            ((n_overline[i] * omega[i] - n_overline[j] * grad_r[j]) * (n_overline[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth(omega) * ds_r + \
            ((n_overline[i] * omega[i] - n_overline[j] * grad_R[j]) * (n_overline[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth(omega) * ds_R
    )
F = F_z + F_omega + F_N

#boundary conditions
bc_r = DirichletBC(Q_z_omega.sub(0), Expression('C', element = Q_z_omega.sub(0).ufl_element(), C = C), boundary_r)
bc_R = DirichletBC(Q_z_omega.sub(0), Expression('2*C', element = Q_z_omega.sub(0).ufl_element(), C = C), boundary_R)
bcs = [bc_r, bc_R]

#solve
solve(F == 0, z_omega, bcs)

#print out the solution
z_dummy, omega_dummy = z_omega.split( deepcopy=True )
xdmffile_z.write( z_dummy, 0 )
xdmffile_omega.write( omega_dummy, 0 )

Here are two meshes: coarse_mesh.zip and fine_mesh.zip, please unzip each of them and put the respective mesh files in [some_path]/coarse_mesh and [some_path]/fine_mesh.
Then if you run on a coarse mesh everything is fine

$ python3 mwe.py /home/fenics/shared/coarse_mesh
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.958e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 6.154e+00 (tol = 1.000e-10) r (rel) = 1.555e-02 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 2.618e-01 (tol = 1.000e-10) r (rel) = 6.615e-04 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 2.722e-04 (tol = 1.000e-10) r (rel) = 6.878e-07 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 1.927e-09 (tol = 1.000e-10) r (rel) = 4.869e-12 (tol = 1.000e-09)
  Newton solver finished in 4 iterations and 4 linear solver iterations.

while if you run on the fine mesh it crashes

$ python3 mwe.py /home/fenics/shared/mesh/membrane_mesh/fine_mesh
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) = 1.366e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

Traceback (most recent call last):
  File "mwe.py", line 189, in <module>
    solve(F == 0, z_omega, bcs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 266, 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 successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

I hope that this helps!

[dokken]

First of all, I would add:

parameters["form_compiler"]["quadrature_degree"] = 5

to not over estimate the quadrature degree of your form.

Secondly, you can set:

solve(F == 0, z_omega, bcs, solver_parameters={
      "newton_solver": {"linear_solver": 'superlu'}})

which will work (slowly):

  Newton iteration 0: r (abs) = 1.366e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
   Newton iteration 1: r (abs) = 1.167e+01 (tol = 1.000e-10) r (rel) = 8.544e-02 (tol = 1.000e-09)
    Newton iteration 2: r (abs) = 1.547e+06 (tol = 1.000e-10) r (rel) = 1.132e+04 (tol = 1.000e-09)

allthough in your case, the solution diverges.
You might have to:

1. Set an initial condition, as Newton is not guaranteed to converge.
2. Use a higher quadrature degree (as underestimating quadrature is a variational crime: https://people.maths.ox.ac.uk/suli/fem.pdf section 3.4).

Choosing degree=10 yields:

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.366e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 5.749e+00 (tol = 1.000e-10) r (rel) = 4.207e-02 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 2.479e-01 (tol = 1.000e-10) r (rel) = 1.814e-03 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 9.497e-04 (tol = 1.000e-10) r (rel) = 6.950e-06 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 1.847e-08 (tol = 1.000e-10) r (rel) = 1.352e-10 (tol = 1.000e-09)
  Newton solver finished in 4 iterations and 4 linear solver iterations.

[mekong]

@dokken

Choosing degree=10 yields

Are you sure you ran it with the fine mesh? I ask because your modifications doe not work for me with the fine_mesh: here is the code where I set the quadrature_degree to 10 and used super_lu as you suggested:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import argparse
import ufl as ufl

dolfin.parameters["form_compiler"]["quadrature_degree"] = 10

####
i, j, k, l = ufl.indices(4)

def my_norm(x):
    return (sqrt(np.dot(x, x)))

def ufl_norm(x):
    return(sqrt(ufl.dot(x, x)))

def X(z):
    x = ufl.SpatialCoordinate(mesh)
    return as_tensor([x[0], x[1], z])

def e(omega):
    return as_tensor([[1, 0, omega[0]], [0, 1, omega[1]]])

def normal(omega):
    return as_tensor(cross(e(omega)[0], e(omega)[1]) /  ufl_norm(cross(e(omega)[0], e(omega)[1])) )

def b(omega):
    return as_tensor((normal(omega))[k] * (e(omega)[i, k]).dx(j), (i,j))

def g(omega):
    return as_tensor([[1+ (omega[0])**2, (omega[0])*(omega[1])],[(omega[0])*(omega[1]), 1+ (omega[1])**2]])

def g_c(omega):
    return ufl.inv(g(omega))

def detg(omega):
    return ufl.det(g(omega))

def abs_detg(omega):
    return np.abs(ufl.det(g(omega)))

def sqrt_detg(omega):
    return sqrt(detg(omega))

def sqrt_abs_detg(omega):
    return sqrt(abs_detg(omega))

def w():
    x = ufl.SpatialCoordinate(mesh)
    return as_tensor([-x[1], x[0]])

def sqrt_deth(omega):
    return(sqrt((w())[i]*(w())[j]*g(omega)[i, j]))

def calc_normal_cg2(mesh):
    n = FacetNormal(mesh)
    V = VectorFunctionSpace(mesh, "CG", 2)
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(u, v) * ds
    l = inner(n, v) * ds
    A = assemble(a, keep_diagonal=True)
    L = assemble(l)

    A.ident_zeros()
    nh = Function(V)
    solve(A, nh.vector(), L)
    return nh

def n(omega):
    u = n_overline
    return as_tensor(u[k]/sqrt(g(omega)[i,j]*u[i]*u[j]), (k))

def H(omega):
    return (0.5 * g_c(omega)[i, j]*b(omega)[j, i])

def K(omega):
    return(ufl.det(as_tensor(b(omega)[i,k]*g_c(omega)[k,j], (i, j))))

def Gamma(omega):
    return as_tensor(0.5 * g_c(omega)[i,l] * ( (g(omega)[l, k]).dx(j) + (g(omega)[j, l]).dx(k) - (g(omega)[j, k]).dx(l) ), (i, j, k))

def Nabla_v(u, omega):
    return as_tensor((u[i]).dx(j) + u[k] * Gamma(omega)[i, k, j], (i, j))

def Nabla_f(f, omega):
    return as_tensor((f[i]).dx(j) - f[k] * Gamma(omega)[k, i, j], (i, j))
####


parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
args = parser.parse_args()

#read mesh
mesh=Mesh()
with XDMFFile((args.input_directory) + "/triangle_mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile((args.input_directory) + "/line_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")

n_overline = FacetNormal( mesh )

# Define boundaries and obstacle
boundary = 'on_boundary'
boundary_r = 'on_boundary && sqrt(pow(x[0], 2) + pow(x[1], 2)) < (1.0+2.0)/2.0'
boundary_R = 'on_boundary && sqrt(pow(x[0], 2) + pow(x[1], 2)) > (1.0+2.0)/2.0'


r = 1.0
R = 2.0
kappa = 1.0
alpha = 1000.0
sigma0 = 1.0
C = 0.5

set_log_level(20)


# Define function spaces
P_z = FiniteElement('P', triangle, 1)
P_omega = VectorElement('P', triangle, 3)
element = MixedElement([P_z, P_omega])
Q_z_omega = FunctionSpace(mesh, element)
Q_z = Q_z_omega.sub(0).collapse()
Q_omega = Q_z_omega.sub(1).collapse()

# Create XDMF files for visualization output
xdmffile_z = XDMFFile('z.xdmf')
xdmffile_omega = XDMFFile('omega.xdmf')

#read an object with label subdomain_id from xdmf file and assign to it the ds `ds_inner`
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_r = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_R = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)

# Define trial and test functions
nu_z, nu_omega = TestFunctions(Q_z_omega)
z_omega = Function(Q_z_omega)
d_z_omega = TrialFunction(Q_z_omega)
sigma = Function(Q_z)
z, omega = split(z_omega)

class grad_r_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = C * x[0] / my_norm( x )
        values[1] = C * x[1] / my_norm( x )
        return (2,)

class grad_R_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = 0
        values[1] = 0
    def value_shape(self):
        return (2,)

class sigma_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = sigma0
    def value_shape(self):
        return (1,)


# Define expressions used in variational forms
kappa = Constant(kappa)
sigma.interpolate(sigma_Expression(element=Q_z.ufl_element()))
grad_r = interpolate(grad_r_Expression(element=Q_omega.ufl_element()), Q_omega)
grad_R = interpolate(grad_R_Expression(element=Q_omega.ufl_element()), Q_omega)

# Define functional for variational problem
F_z = ( kappa * ( g_c(omega)[i, j] * (H(omega).dx(j)) * (nu_z.dx(i)) - 2.0 * H(omega) * ( (H(omega))**2 - K(omega) ) * nu_z ) + sigma * H(omega) * nu_z ) * sqrt_detg(omega) * dx \
    - ( kappa * (n(omega))[i] * nu_z * (H(omega).dx(i)) ) * sqrt_deth(omega) * ds
F_omega = ( - z * Nabla_v(nu_omega, omega)[i, i] - omega[i] * nu_omega[i] ) *  sqrt_detg(omega) * dx + \
          ( (n(omega))[i] * g(omega)[i, j] * z * nu_omega[j] ) * sqrt_deth(omega) * ds
#Nitche's part of the functional
F_N = alpha * ( \
            ((n_overline[i] * omega[i] - n_overline[j] * grad_r[j]) * (n_overline[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth(omega) * ds_r + \
            ((n_overline[i] * omega[i] - n_overline[j] * grad_R[j]) * (n_overline[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth(omega) * ds_R
    )
F = F_z + F_omega + F_N

#boundary conditions
bc_r = DirichletBC(Q_z_omega.sub(0), Expression('C', element = Q_z_omega.sub(0).ufl_element(), C = C), boundary_r)
bc_R = DirichletBC(Q_z_omega.sub(0), Expression('2*C', element = Q_z_omega.sub(0).ufl_element(), C = C), boundary_R)
bcs = [bc_r, bc_R]

#solve
J  = derivative(F, z_omega, d_z_omega)  # Gateaux derivative in dir. of du
problem = NonlinearVariationalProblem(F, z_omega, bcs, J)
solver  = NonlinearVariationalSolver(problem)
params = {'nonlinear_solver': 'newton', 'newton_solver': { 'linear_solver' : 'superlu' }, }
solver.parameters.update(params)
solver.solve()

#print out the solution
z_dummy, omega_dummy = z_omega.split( deepcopy=True )
xdmffile_z.write( z_dummy, 0 )
xdmffile_omega.write( omega_dummy, 0 )

it runs fine with the coarse mesh

$  python3 mwe.py /home/fenics/shared/mesh/membrane_mesh/coarse_mesh
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.958e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 6.152e+00 (tol = 1.000e-10) r (rel) = 1.555e-02 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 2.609e-01 (tol = 1.000e-10) r (rel) = 6.593e-04 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 2.790e-04 (tol = 1.000e-10) r (rel) = 7.049e-07 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 2.050e-09 (tol = 1.000e-10) r (rel) = 5.180e-12 (tol = 1.000e-09)
  Newton solver finished in 4 iterations and 4 linear solver iterations.

but it crashes with the fine mesh

$python3 mwe.py /home/fenics/shared/mesh/membrane_mesh/fine_mesh
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.366e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Killed

Any idea ?

[dokken]

I ran on the fine mesh, as the coarse mesh link didn’t work for me.

It seems like you are running out of memory on your system, as this is what kill indicates.

Try running in parallel with superlu_dist instead of superlu

[mekong]

If I replace superlu with superlu_dist as follows:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import argparse
import ufl as ufl

dolfin.parameters["form_compiler"]["quadrature_degree"] = 10

####
i, j, k, l = ufl.indices(4)

def my_norm(x):
    return (sqrt(np.dot(x, x)))

def ufl_norm(x):
    return(sqrt(ufl.dot(x, x)))

def X(z):
    x = ufl.SpatialCoordinate(mesh)
    return as_tensor([x[0], x[1], z])

def e(omega):
    return as_tensor([[1, 0, omega[0]], [0, 1, omega[1]]])

def normal(omega):
    return as_tensor(cross(e(omega)[0], e(omega)[1]) /  ufl_norm(cross(e(omega)[0], e(omega)[1])) )

def b(omega):
    return as_tensor((normal(omega))[k] * (e(omega)[i, k]).dx(j), (i,j))

def g(omega):
    return as_tensor([[1+ (omega[0])**2, (omega[0])*(omega[1])],[(omega[0])*(omega[1]), 1+ (omega[1])**2]])

def g_c(omega):
    return ufl.inv(g(omega))

def detg(omega):
    return ufl.det(g(omega))

def abs_detg(omega):
    return np.abs(ufl.det(g(omega)))

def sqrt_detg(omega):
    return sqrt(detg(omega))

def sqrt_abs_detg(omega):
    return sqrt(abs_detg(omega))

def w():
    x = ufl.SpatialCoordinate(mesh)
    return as_tensor([-x[1], x[0]])

def sqrt_deth(omega):
    return(sqrt((w())[i]*(w())[j]*g(omega)[i, j]))

def calc_normal_cg2(mesh):
    n = FacetNormal(mesh)
    V = VectorFunctionSpace(mesh, "CG", 2)
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(u, v) * ds
    l = inner(n, v) * ds
    A = assemble(a, keep_diagonal=True)
    L = assemble(l)

    A.ident_zeros()
    nh = Function(V)
    solve(A, nh.vector(), L)
    return nh

def n(omega):
    u = n_overline
    return as_tensor(u[k]/sqrt(g(omega)[i,j]*u[i]*u[j]), (k))

def H(omega):
    return (0.5 * g_c(omega)[i, j]*b(omega)[j, i])

def K(omega):
    return(ufl.det(as_tensor(b(omega)[i,k]*g_c(omega)[k,j], (i, j))))

def Gamma(omega):
    return as_tensor(0.5 * g_c(omega)[i,l] * ( (g(omega)[l, k]).dx(j) + (g(omega)[j, l]).dx(k) - (g(omega)[j, k]).dx(l) ), (i, j, k))

def Nabla_v(u, omega):
    return as_tensor((u[i]).dx(j) + u[k] * Gamma(omega)[i, k, j], (i, j))

def Nabla_f(f, omega):
    return as_tensor((f[i]).dx(j) - f[k] * Gamma(omega)[k, i, j], (i, j))
####


parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
args = parser.parse_args()

#read mesh
mesh=Mesh()
with XDMFFile((args.input_directory) + "/triangle_mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile((args.input_directory) + "/line_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")

n_overline = FacetNormal( mesh )

# Define boundaries and obstacle
boundary = 'on_boundary'
boundary_r = 'on_boundary && sqrt(pow(x[0], 2) + pow(x[1], 2)) < (1.0+2.0)/2.0'
boundary_R = 'on_boundary && sqrt(pow(x[0], 2) + pow(x[1], 2)) > (1.0+2.0)/2.0'


r = 1.0
R = 2.0
kappa = 1.0
alpha = 1000.0
sigma0 = 1.0
C = 0.5

set_log_level(20)


# Define function spaces
P_z = FiniteElement('P', triangle, 1)
P_omega = VectorElement('P', triangle, 3)
element = MixedElement([P_z, P_omega])
Q_z_omega = FunctionSpace(mesh, element)
Q_z = Q_z_omega.sub(0).collapse()
Q_omega = Q_z_omega.sub(1).collapse()

# Create XDMF files for visualization output
xdmffile_z = XDMFFile('z.xdmf')
xdmffile_omega = XDMFFile('omega.xdmf')

#read an object with label subdomain_id from xdmf file and assign to it the ds `ds_inner`
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_r = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_R = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)

# Define trial and test functions
nu_z, nu_omega = TestFunctions(Q_z_omega)
z_omega = Function(Q_z_omega)
d_z_omega = TrialFunction(Q_z_omega)
sigma = Function(Q_z)
z, omega = split(z_omega)

class grad_r_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = C * x[0] / my_norm( x )
        values[1] = C * x[1] / my_norm( x )
        return (2,)

class grad_R_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = 0
        values[1] = 0
    def value_shape(self):
        return (2,)

class sigma_Expression( UserExpression ):
    def eval(self, values, x):
        values[0] = sigma0
    def value_shape(self):
        return (1,)


# Define expressions used in variational forms
kappa = Constant(kappa)
sigma.interpolate(sigma_Expression(element=Q_z.ufl_element()))
grad_r = interpolate(grad_r_Expression(element=Q_omega.ufl_element()), Q_omega)
grad_R = interpolate(grad_R_Expression(element=Q_omega.ufl_element()), Q_omega)

# Define functional for variational problem
F_z = ( kappa * ( g_c(omega)[i, j] * (H(omega).dx(j)) * (nu_z.dx(i)) - 2.0 * H(omega) * ( (H(omega))**2 - K(omega) ) * nu_z ) + sigma * H(omega) * nu_z ) * sqrt_detg(omega) * dx \
    - ( kappa * (n(omega))[i] * nu_z * (H(omega).dx(i)) ) * sqrt_deth(omega) * ds
F_omega = ( - z * Nabla_v(nu_omega, omega)[i, i] - omega[i] * nu_omega[i] ) *  sqrt_detg(omega) * dx + \
          ( (n(omega))[i] * g(omega)[i, j] * z * nu_omega[j] ) * sqrt_deth(omega) * ds
#Nitche's part of the functional
F_N = alpha * ( \
            ((n_overline[i] * omega[i] - n_overline[j] * grad_r[j]) * (n_overline[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth(omega) * ds_r + \
            ((n_overline[i] * omega[i] - n_overline[j] * grad_R[j]) * (n_overline[k] * g( omega )[k, l] * nu_omega[l])) * sqrt_deth(omega) * ds_R
    )
F = F_z + F_omega + F_N

#boundary conditions
bc_r = DirichletBC(Q_z_omega.sub(0), Expression('C', element = Q_z_omega.sub(0).ufl_element(), C = C), boundary_r)
bc_R = DirichletBC(Q_z_omega.sub(0), Expression('2*C', element = Q_z_omega.sub(0).ufl_element(), C = C), boundary_R)
bcs = [bc_r, bc_R]

#solve
J  = derivative(F, z_omega, d_z_omega)  # Gateaux derivative in dir. of du
problem = NonlinearVariationalProblem(F, z_omega, bcs, J)
solver  = NonlinearVariationalSolver(problem)
params = {'nonlinear_solver': 'newton', 'newton_solver': { 'linear_solver' : 'superlu_dist' }, }
solver.parameters.update(params)
solver.solve()

#print out the solution
z_dummy, omega_dummy = z_omega.split( deepcopy=True )
xdmffile_z.write( z_dummy, 0 )
xdmffile_omega.write( omega_dummy, 0 )

and run, I get

mpirun -np 6 python3 mwe.py /home/fenics/shared/mesh/membrane_mesh/fine_mesh

Process 0: Solving nonlinear variational problem.
Process 2: Solving nonlinear variational problem.Process 4: Solving nonlinear variational problem.

Process 5: Solving nonlinear variational problem.Process 3: Solving nonlinear variational problem.

Process 1: Solving nonlinear variational problem.
Traceback (most recent call last):
  File "mwe.py", line 201, 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 linear system.
*** Reason:  Unknown solver method "superlu_dist". Use list_linear_solver_methods() to list available methods.
*** Where:   This error was encountered inside LinearSolver.cpp.
*** Process: 4
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Traceback (most recent call last):
  File "mwe.py", line 201, 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 linear system.
*** Reason:  Unknown solver method "superlu_dist". Use list_linear_solver_methods() to list available methods.
*** Where:   This error was encountered inside LinearSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Traceback (most recent call last):
  File "mwe.py", line 201, 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 linear system.
*** Reason:  Unknown solver method "superlu_dist". Use list_linear_solver_methods() to list available methods.
*** Where:   This error was encountered inside LinearSolver.cpp.
*** Process: 5
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Traceback (most recent call last):
  File "mwe.py", line 201, 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 linear system.
*** Reason:  Unknown solver method "superlu_dist". Use list_linear_solver_methods() to list available methods.
*** Where:   This error was encountered inside LinearSolver.cpp.
*** Process: 1
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Traceback (most recent call last):
  File "mwe.py", line 201, 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 linear system.
*** Reason:  Unknown solver method "superlu_dist". Use list_linear_solver_methods() to list available methods.
*** Where:   This error was encountered inside LinearSolver.cpp.
*** Process: 2
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Traceback (most recent call last):
  File "mwe.py", line 201, 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 linear system.
*** Reason:  Unknown solver method "superlu_dist". Use list_linear_solver_methods() to list available methods.
*** Where:   This error was encountered inside LinearSolver.cpp.
*** Process: 3
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

[mekong]

Hello, is there any update about my last post on this topic?

Thank you

[dokken]

Well, it doesn’t seem like your installation of legacy fenics has superlu_dist installed.

As previously stated, it seems like your computer doesn’t have sufficient memory to run the finer simulation.

[mekong]

How may I install fenics with superlu_dist ? (Yes, I did look for this topic on the forum, none of the solutions work for me, see this or this).

[Stein]

Sorry that this is going in a completely different direction, but have you tried using an iterative solver? When problems get too large and too expensive for direct solvers (mumps, superlu), the natural choice seems to move to iterative solvers (CG/ByCGstab/GMRes). Preconditioning may be a bit tricky for your mixed problem though.

[mekong]

@Stein

the natural choice seems to move to iterative solvers (CG/ByCGstab/GMRes). Preconditioning may be a bit tricky for your mixed problem though.

Thankyou. I tried ‘cg’ and ‘gmres’ and they both give the error

*** 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 linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 2 iterations (PETSc reason DIVERGED_INDEFINITE_PC, residual norm ||r|| = 4.666663e+40).

*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f

I tried ‘bycgstab’ but it is not recognized as a solver.