Error in Linear Variational Solver

When trying to run the program “ft10_poisson_extended.py” which can be found in

https://raw.githubusercontent.com/hplgit/fenics-tutorial/master/pub/python/vol1/ft10_poisson_extended.py

I received this error

prm = LinearVariationalSolver.default_parameters()

AttributeError: type object ‘dolfin.cpp.fem.LinearVariationalSolver’ has no attribute ‘default_parameters’

How to run this program in order to get started with it?

Basically I am trying to run the demos in page. 98 of tutorial, vol. 1
from the tutorial in https://fenicsproject.org/tutorial/

Note that the fenics tutorial was published in 2017, to be compatible with the 2016.2.0 version. To be able to help you, we need to know which version of FEniCS you are using.
For your particular issue, see for instance: Updating Fenics Vol 1 Demos

@dokken:

  1. I am using Spyder (Python 3.7) . Although I cannot run the command in Spyder, but in Linux terminal, I ran

dolfin-version

I received

2018.1.0

  1. I replaced the following

    prm = LinearVariationalSolver.default_parameters()
    if linear_solver == ‘Krylov’:
    prm.linear_solver = ‘gmres’
    prm.preconditioner = ‘ilu’
    prm.krylov_solver.absolute_tolerance = abs_tol
    prm.krylov_solver.relative_tolerance = rel_tol
    prm.krylov_solver.maximum_iterations = max_iter
    else:
    prm.linear_solver = ‘lu’

with your recommendation

prm = LinearVariationalSolver.default_parameters()
if linear_solver == 'Krylov':
    prm["linear_solver"] = 'gmres'
    prm["preconditioner"] = 'ilu'
    prm["krylov_solver.absolute_tolerance"] = abs_tol
    prm["krylov_solver.relative_tolerance"] = rel_tol
    prm["krylov_solver.maximum_iterations"] = max_iter
else:
    prm["linear_solver"] = 'lu'

But still the receive the error

demos[nr]()

TypeError: tuple indices must be integers or slices, not str

After runnign the python program, I use
demos[0]()

How to make it run?

As this demo is created for python2.7, and especially dolfin 2016.2.0, I would use that version to run the tutorials. If you would like to use a later version of dolfin for the tutorials, you will have to convert the tutorial to python3 syntax. Most of the changes from python2.7 to 3 can be done automatically by using:
2to3 filename.py -w.
However this demo (with the specific parameters) will not run with dolfin 2018.1.0.
You should use the latest stable release:

2019.1.0

Note that

is not the correct syntax, and should be:

 prm["krylov_solver"]["absolute_tolerance"] = abs_tol

Here is a converted version of the tutorial that runs with the latest release of dolfin.

"""
FEniCS tutorial demo program: Poisson equation with a combination of
Dirichlet, Neumann, and Robin boundary conditions.

  -div(kappa*grad(u)) = f

This program illustrates a number of different topics:

- How to solve a problem using three different approaches of varying
  complexity: solve / LinearVariationalSolver / assemble + solve
- How to compute fluxes
- How to set combinations of boundary conditions
- How to set parameters for linear solvers
- How to create manufactured solutions with SymPy
- How to create unit tests
- How to represent solutions as structured fields
"""



from fenics import *
from boxfield import *
import numpy as np

#---------------------------------------------------------------------
# Solvers
#---------------------------------------------------------------------

def solver(kappa, f, u_D, Nx, Ny,
           degree=1,
           linear_solver='Krylov',
           abs_tol=1E-5,
           rel_tol=1E-3,
           max_iter=1000):
    """
    Solve -div(kappa*grad(u) = f on (0, 1) x (0, 1) with 2*Nx*Ny Lagrange
    elements of specified degree and u = u_D on the boundary.
    """

    # Create mesh and define function space
    mesh = UnitSquareMesh(Nx, Ny)
    V = FunctionSpace(mesh, 'P', degree)

    # Define boundary condition
    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = kappa*dot(grad(u), grad(v))*dx
    L = f*v*dx

    # Set linear solver parameters
    prm = LinearVariationalSolver.default_parameters()
    if linear_solver == 'Krylov':
        prm["linear_solver"] = 'gmres'
        prm["preconditioner"] = 'ilu'
        prm["krylov_solver"]["absolute_tolerance"] = abs_tol
        prm["krylov_solver"]["relative_tolerance"] = rel_tol
        prm["krylov_solver"]["maximum_iterations"] = max_iter
    else:
        prm["linear_solver"] = 'lu'

    # Compute solution
    u = Function(V)
    solve(a == L, u, bc, solver_parameters=prm)

    return u

def solver_objects(kappa, f, u_D, Nx, Ny,
                   degree=1,
                   linear_solver='Krylov',
                   abs_tol=1E-5,
                   rel_tol=1E-3,
                   max_iter=1000):
    "Same as the solver() function but using LinearVariationalSolver"

    # Create mesh and define function space
    mesh = UnitSquareMesh(Nx, Ny)
    V = FunctionSpace(mesh, 'P', degree)

    # Define boundary condition
    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = kappa*dot(grad(u), grad(v))*dx
    L = f*v*dx

    # Compute solution
    u = Function(V)
    problem = LinearVariationalProblem(a, L, u, bc)
    solver = LinearVariationalSolver(problem)

    # Set linear solver parameters
    prm = solver.parameters
    if linear_solver == 'Krylov':
        prm.linear_solver = 'gmres'
        prm.preconditioner = 'ilu'
        prm.krylov_solver.absolute_tolerance = abs_tol
        prm.krylov_solver.relative_tolerance = rel_tol
        prm.krylov_solver.maximum_iterations = max_iter
    else:
        prm.linear_solver = 'lu'

    # Compute solution
    solver.solve()

    return u

def solver_linalg(kappa, f, u_D, Nx, Ny,
                 degree=1,
                 linear_solver='Krylov',
                 abs_tol=1E-5,
                 rel_tol=1E-3,
                 max_iter=1000):
    "Same as the solver() function but assembling and solving Ax = b"

    # Create mesh and define function space
    mesh = UnitSquareMesh(Nx, Ny)
    V = FunctionSpace(mesh, 'P', degree)

    # Define boundary condition
    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = kappa*dot(grad(u), grad(v))*dx
    L = f*v*dx

    # Assemble linear system
    A = assemble(a)
    b = assemble(L)

    # Apply boundary conditions
    bc.apply(A, b)

    # Create linear solver and set parameters
    if linear_solver == 'Krylov':
        solver = KrylovSolver('gmres', 'ilu')
        solver.parameters.absolute_tolerance = abs_tol
        solver.parameters.relative_tolerance = rel_tol
        solver.parameters.maximum_iterations = max_iter
    else:
        solver = LUSolver()

    # Compute solution
    u = Function(V)
    solver.solve(A, u.vector(), b)

    return u

def solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
               degree=1,
               subdomains=[],
               linear_solver='Krylov',
               abs_tol=1E-5,
               rel_tol=1E-3,
               max_iter=1000):
    """
    Solve -div(kappa*grad(u) = f on (0, 1) x (0, 1) with 2*Nx*Ny Lagrange
    elements of specified degree and u = u_D on the boundary. This version
    of the solver uses a specified combination of Dirichlet, Neumann, and
    Robin boundary conditions.
    """

    # Create mesh and define function space
    mesh = UnitSquareMesh(Nx, Ny)
    V = FunctionSpace(mesh, 'P', degree)

    # Check if we have subdomains
    if subdomains:
        if not isinstance(kappa, (list, tuple, np.ndarray)):
            raise TypeError(
                'kappa must be array if we have sudomains, not %s'
                % type(kappa))
        materials = CellFunction('size_t', mesh)
        materials.set_all(0)
        for m, subdomain in enumerate(subdomains[1:], 1):
            subdomain.mark(materials, m)

        kappa_values = kappa
        V0 = FunctionSpace(mesh, 'DG', 0)
        kappa  = Function(V0)
        help = np.asarray(materials.array(), dtype=np.int32)
        kappa.vector()[:] = np.choose(help, kappa_values)
    else:
        if not isinstance(kappa, (Expression, Constant)):
            raise TypeError(
                'kappa is type %s, must be Expression or Constant'
                % type(kappa))

    # Define boundary subdomains
    tol = 1e-14

    class BoundaryX0(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0], 0, tol)

    class BoundaryX1(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0], 1, tol)

    class BoundaryY0(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[1], 0, tol)

    class BoundaryY1(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[1], 1, tol)

    # Mark boundaries
    boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
    boundary_markers.set_all(9999)
    bx0 = BoundaryX0()
    bx1 = BoundaryX1()
    by0 = BoundaryY0()
    by1 = BoundaryY1()
    bx0.mark(boundary_markers, 0)
    bx1.mark(boundary_markers, 1)
    by0.mark(boundary_markers, 2)
    by1.mark(boundary_markers, 3)

    # Redefine boundary integration measure
    ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

    # Collect Dirichlet conditions
    bcs = []
    for i in boundary_conditions:
        if 'Dirichlet' in boundary_conditions[i]:
            bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                             boundary_markers, i)
            bcs.append(bc)

    debug1 = False
    if debug1:

        # Print all vertices that belong to the boundary parts
        for x in mesh.coordinates():
            if bx0.inside(x, True): print('%s is on x = 0' % x)
            if bx1.inside(x, True): print('%s is on x = 1' % x)
            if by0.inside(x, True): print('%s is on y = 0' % x)
            if by1.inside(x, True): print('%s is on y = 1' % x)

        # Print the Dirichlet conditions
        print('Number of Dirichlet conditions:', len(bcs))
        if V.ufl_element().degree() == 1:  # P1 elements
            d2v = dof_to_vertex_map(V)
            coor = mesh.coordinates()
            for i, bc in enumerate(bcs):
                print('Dirichlet condition %d' % i)
                boundary_values = bc.get_boundary_values()
                for dof in boundary_values:
                    print('   dof %2d: u = %g' % (dof, boundary_values[dof]))
                    if V.ufl_element().degree() == 1:
                        print('    at point %s' %
                              (str(tuple(coor[d2v[dof]].tolist()))))

    # Define trial and test functions
    u = TrialFunction(V)
    v = TestFunction(V)

    # Collect Neumann integrals
    integrals_N = []
    for i in boundary_conditions:
        if 'Neumann' in boundary_conditions[i]:
            if boundary_conditions[i]['Neumann'] != 0:
                g = boundary_conditions[i]['Neumann']
                integrals_N.append(g*v*ds(i))

    # Collect Robin integrals
    integrals_R_a = []
    integrals_R_L = []
    for i in boundary_conditions:
        if 'Robin' in boundary_conditions[i]:
            r, s = boundary_conditions[i]['Robin']
            integrals_R_a.append(r*u*v*ds(i))
            integrals_R_L.append(r*s*v*ds(i))

    # Simpler Robin integrals
    integrals_R = []
    for i in boundary_conditions:
        if 'Robin' in boundary_conditions[i]:
            r, s = boundary_conditions[i]['Robin']
            integrals_R.append(r*(u - s)*v*ds(i))

    # Sum integrals to define variational problem
    a = kappa*dot(grad(u), grad(v))*dx + sum(integrals_R_a)
    L = f*v*dx - sum(integrals_N) + sum(integrals_R_L)

    # Simpler variational problem
    F = kappa*dot(grad(u), grad(v))*dx + \
        sum(integrals_R) - f*v*dx + sum(integrals_N)
    a, L = lhs(F), rhs(F)

    # Set linear solver parameters
    prm = LinearVariationalSolver.default_parameters()
    if linear_solver == 'Krylov':
        prm["linear_solver"] = 'gmres'
        prm["preconditioner"] = 'ilu'
        prm["krylov_solver"]["absolute_tolerance"] = abs_tol
        prm["krylov_solver"]["relative_tolerance"] = rel_tol
        prm["krylov_solver"]["maximum_iterations"] = max_iter
    else:
        prm["linear_solver"] = 'lu'

    # Compute solution
    u = Function(V)
    solve(a == L, u, bcs, solver_parameters=prm)

    return u

#---------------------------------------------------------------------
# Utility functions
#---------------------------------------------------------------------

def compute_errors(u_e, u):
    """Compute various measures of the error u - u_e, where
    u is a finite element Function and u_e is an Expression."""

    # Get function space
    V = u.function_space()

    # Explicit computation of L2 norm
    error = (u - u_e)**2*dx
    E1 = sqrt(abs(assemble(error)))

    # Explicit interpolation of u_e onto the same space as u
    u_e_ = interpolate(u_e, V)
    error = (u - u_e_)**2*dx
    E2 = sqrt(abs(assemble(error)))

    # Explicit interpolation of u_e to higher-order elements.
    # u will also be interpolated to the space Ve before integration
    Ve = FunctionSpace(V.mesh(), 'P', 5)
    u_e_ = interpolate(u_e, Ve)
    error = (u - u_e)**2*dx
    E3 = sqrt(abs(assemble(error)))

    # Infinity norm based on nodal values
    u_e_ = interpolate(u_e, V)
    E4 = abs(u_e_.vector().get_local() - u.vector().get_local()).max()

    # L2 norm
    E5 = errornorm(u_e, u, norm_type='L2', degree_rise=3)

    # H1 seminorm
    E6 = errornorm(u_e, u, norm_type='H10', degree_rise=3)

    # Collect error measures in a dictionary with self-explanatory keys
    errors = {'u - u_e': E1,
              'u - interpolate(u_e, V)': E2,
              'interpolate(u, Ve) - interpolate(u_e, Ve)': E3,
              'infinity norm (of dofs)': E4,
              'L2 norm': E5,
              'H10 seminorm': E6}

    return errors

def compute_convergence_rates(u_e, f, u_D, kappa,
                              max_degree=3, num_levels=5):
    "Compute convergences rates for various error norms"

    h = {}  # discretization parameter: h[degree][level]
    E = {}  # error measure(s): E[degree][level][error_type]

    # Iterate over degrees and mesh refinement levels
    degrees = list(range(1, max_degree + 1))
    for degree in degrees:
        n = 8  # coarsest mesh division
        h[degree] = []
        E[degree] = []
        for i in range(num_levels):
            h[degree].append(1.0 / n)
            u = solver(kappa, f, u_D, n, n, degree, linear_solver='direct')
            errors = compute_errors(u_e, u)
            E[degree].append(errors)
            print('2 x (%d x %d) P%d mesh, %d unknowns, E1 = %g' %
              (n, n, degree, u.function_space().dim(), errors['u - u_e']))
            n *= 2

    # Compute convergence rates
    from math import log as ln  # log is a fenics name too
    etypes = list(E[1][0].keys())
    rates = {}
    for degree in degrees:
        rates[degree] = {}
        for error_type in sorted(etypes):
            rates[degree][error_type] = []
            for i in range(1, num_levels):
                Ei = E[degree][i][error_type]
                Eim1 = E[degree][i - 1][error_type]
                r = ln(Ei / Eim1) / ln(h[degree][i] / h[degree][i - 1])
                rates[degree][error_type].append(round(r, 2))

    return etypes, degrees, rates

def flux(u, kappa):
    "Return -kappa*grad(u) projected into same space as u"
    V = u.function_space()
    mesh = V.mesh()
    degree = V.ufl_element().degree()
    W = VectorFunctionSpace(mesh, 'P', degree)
    flux_u = project(-kappa*grad(u), W)
    return flux_u

def normalize_solution(u):
    "Normalize u: return u divided by max(u)"
    u_array = u.vector().get_local()
    u_max = np.max(np.abs(u_array))
    u_array /= u_max
    u.vector()[:] = u_array
    #u.vector().set_local(u_array)  # alternative
    return u

#---------------------------------------------------------------------
# Unit tests (functions beginning with test_)
# These unit tests can be run by calling `py.test poisson_extended.py`
#---------------------------------------------------------------------

def test_solvers():
    "Reproduce exact solution to machine precision with different solvers"

    solver_functions = (solver, solver_objects, solver_linalg)
    tol = {'direct': {1: 1E-11, 2: 1E-11, 3: 1E-11},
           'Krylov': {1: 1E-14, 2: 1E-05, 3: 1E-03}}
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
    kappa = Expression('x[0] + x[1]', degree=1)
    f = Expression('-8*x[0] - 10*x[1]', degree=1)
    for Nx, Ny in [(3, 3), (3, 5), (5 ,3)]:
        for degree in 1, 2, 3:
            for linear_solver in 'direct', 'Krylov':
                for solver_function in solver_functions:
                    print('solving on 2 x (%d x %d) mesh with P%d elements'
                          % (Nx, Ny, degree)),
                    print(' %s solver, %s function' %
                          (linear_solver, solver_function.__name__))
                    u = solver_function(kappa, f, u_D, Nx, Ny, degree,
                                        linear_solver=linear_solver,
                                        abs_tol=0.1*tol[linear_solver][degree],
                                        rel_tol=0.1*tol[linear_solver][degree])
                    V = u.function_space()
                    u_D_Function = interpolate(u_D, V)
                    u_D_array = u_D_Function.vector().get_local()
                    error_max = (u_D_array - u.vector().get_local()).max()
                    msg = 'max error: %g for 2 x (%d x %d) mesh, ' \
                          'degree = %d, %s solver, %s' % \
                          (error_max, Nx, Ny, degree, linear_solver,
                           solver_function.__name__)
                    print(msg)
                    assert error_max < tol[linear_solver][degree], msg

def test_normalize_solution():
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
    f = Constant(-6.0)
    u = solver(f, u_D, 4, 2, 1, linear_solver='direct')
    u = normalize_solution(u)
    computed = u.vector().get_local().max()
    expected = 1.0
    assert abs(expected - computed) < 1E-15

#---------------------------------------------------------------------
# Demo programs
#---------------------------------------------------------------------

def demo_test():
    "Solve test problem and plot solution"
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
    kappa = Expression('x[0] + x[1]', degree=1)
    f = Expression('-8*x[0] - 10*x[1]', degree=1)
    u = solver(kappa, f, u_D, 8, 8, 1)
    vtkfile = File('poisson_extended/solution_test.pvd')
    vtkfile << u
    plot(u)

def demo_flux(Nx=8, Ny=8):
    "Solve test problem and compute flux"

    # Compute solution
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
    kappa = Expression('x[0] + x[1]', degree=1)
    f = Expression('-8*x[0] - 10*x[1]', degree=1)
    u = solver(kappa, f, u_D, Nx, Ny, 1, linear_solver='direct')

    # Compute and plot flux
    flux_u = flux(u, kappa)
    flux_u_x, flux_u_y = flux_u.split(deepcopy=True)
    plot(u, title=u.name())
    plot(flux_u, title=flux_u.name())
    plot(flux_u_x, title=flux_u_x.name())
    plot(flux_u_y, title=flux_u_y.name())

    # Exact flux expressions
    u_e = lambda x, y: 1 + x**2 + 2*y**2
    flux_x_exact = lambda x, y: -(x+y)*2*x
    flux_y_exact = lambda x, y: -(x+y)*4*y

    # Compute error in flux
    coor = u.function_space().mesh().coordinates()
    for i, value in enumerate(flux_u_x.compute_vertex_values()):
        print('vertex %d, x = %s, -p*u_x = %g, error = %g' %
              (i, tuple(coor[i]), value, flux_x_exact(*coor[i]) - value))
    for i, value in enumerate(flux_u_y.compute_vertex_values()):
        print('vertex %d, x = %s, -p*u_y = %g, error = %g' %
              (i, tuple(coor[i]), value, flux_y_exact(*coor[i]) - value))

def demo_convergence_rates():
    "Compute convergence rates in various norms for P1, P2, P3"

    # Define exact solution and coefficients
    omega = 1.0
    u_e = Expression('sin(omega*pi*x[0])*sin(omega*pi*x[1])',
                     degree=6, omega=omega)
    f = 2*omega**2*pi**2*u_e
    u_D = Constant(0)
    kappa = Constant(1)

    # Compute and print convergence rates
    etypes, degrees, rates = compute_convergence_rates(u_e, f, u_D, kappa)
    for error_type in etypes:
        print('\n' + error_type)
        for degree in degrees:
            print('P%d: %s' % (degree, str(rates[degree][error_type])[1:-1]))

def demo_structured_mesh():
    "Use structured mesh data to create plots with Matplotlib"

    # Define exact solution (Mexican hat) and coefficients
    from sympy import exp, sin, pi
    import sympy as sym
    H = lambda x: exp(-16*(x-0.5)**2)*sin(3*pi*x)
    x, y = sym.symbols('x[0], x[1]')
    u = H(x)*H(y)
    u_code = sym.printing.ccode(u)
    u_code = u_code.replace('M_PI', 'pi')
    print('C code for u:', u_code)
    u_D = Expression(u_code, degree=1)
    kappa = 1  # Note: Can't use Constant(1) here because of sym.diff (!)
    f = sym.diff(-kappa*sym.diff(u, x), x) + \
        sym.diff(-kappa*sym.diff(u, y), y)
    f = sym.simplify(f)
    f_code = sym.printing.ccode(f)
    f_code = f_code.replace('M_PI', 'pi')
    f = Expression(f_code, degree=1)
    flux_u_x_exact = sym.lambdify([x, y], -kappa*sym.diff(u, x),
                                  modules='numpy')
    print('C code for f:', f_code)
    kappa = Constant(1)
    nx = 22;  ny = 22

    # Compute solution and represent as a structured field
    u = solver(kappa, f, u_D, nx, ny, 1, linear_solver='direct')
    u_box = FEniCSBoxField(u, (nx, ny))

    # Set coordinates and extract values
    X = 0; Y = 1
    u_ = u_box.values

    # Iterate over 2D mesh points (i, j)
    debug2 = False
    if debug2:
        for j in range(u_.shape[1]):
            for i in range(u_.shape[0]):
                print('u[%d, %d] = u(%g, %g) = %g' %
                      (i, j,
                       u_box.grid.coor[X][i], u_box.grid.coor[Y][j],
                       u_[i, j]))

    # Make surface plot
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cm
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    cv = u_box.grid.coorv
    ax.plot_surface(cv[X], cv[Y], u_, cmap=cm.coolwarm,
                    rstride=1, cstride=1)
    plt.title('Surface plot of solution')
    plt.savefig('poisson_extended/surface_plot.png')
    plt.savefig('poisson_extended/surface_plot.pdf')

    # Make contour plot
    fig = plt.figure()
    ax = fig.gca()
    cs = ax.contour(cv[X], cv[Y], u_, 7)
    plt.clabel(cs)
    plt.axis('equal')
    plt.title('Contour plot of solution')
    plt.savefig('poisson_extended/contour_plot.png')
    plt.savefig('poisson_extended/contour_plot.pdf')

    # Plot u along a line y = const and compare with exact solution
    start = (0, 0.4)
    x, u_val, y_fixed, snapped = u_box.gridline(start, direction=X)
    u_e_val = [u_D((x_, y_fixed)) for x_ in x]
    plt.figure()
    plt.plot(x, u_val, 'r-')
    plt.plot(x, u_e_val, 'bo')
    plt.legend(['P1 elements', 'exact'], loc='best')
    plt.title('Solution along line y=%g' % y_fixed)
    plt.xlabel('x');  plt.ylabel('u')
    plt.savefig('poisson_extended/line_plot.png')
    plt.savefig('poisson_extended/line_plot.pdf')

    # Plot the numerical and exact flux along the same line
    flux_u = flux(u, kappa)
    flux_u_x, flux_u_y = flux_u.split(deepcopy=True)
    flux2_x = flux_u_x if flux_u_x.ufl_element().degree() == 1 \
              else interpolate(flux_x,
                   FunctionSpace(u.function_space().mesh(), 'P', 1))
    flux_u_x_box = FEniCSBoxField(flux_u_x, (nx,ny))
    x, flux_u_val, y_fixed, snapped = \
       flux_u_x_box.gridline(start, direction=X)
    y = y_fixed
    plt.figure()
    plt.plot(x, flux_u_val, 'r-')
    plt.plot(x, flux_u_x_exact(x, y_fixed), 'bo')
    plt.legend(['P1 elements', 'exact'], loc='best')
    plt.title('Flux along line y=%g' % y_fixed)
    plt.xlabel('x');  plt.ylabel('u')
    plt.savefig('poisson_extended/line_flux.png')
    plt.savefig('poisson_extended/line_flux.pdf')

    plt.show()

def demo_bcs():
    "Compute and plot solution using a combination of boundary conditions"

    # Define manufactured solution in sympy and derive f, g, etc.
    import sympy as sym
    x, y = sym.symbols('x[0], x[1]')            # needed by UFL
    u = 1 + x**2 + 2*y**2                       # exact solution
    u_e = u                                     # exact solution
    u_00 = u.subs(x, 0)                         # restrict to x = 0
    u_01 = u.subs(x, 1)                         # restrict to x = 1
    f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)  # -Laplace(u)
    f = sym.simplify(f)                         # simplify f
    g = -sym.diff(u, y).subs(y, 1)              # compute g = -du/dn
    r = 1000                                    # Robin data, arbitrary
    s = u                                       # Robin data, u = s

    # Collect variables
    variables = [u_e, u_00, u_01, f, g, r, s]

    # Turn into C/C++ code strings
    variables = [sym.printing.ccode(var) for var in variables]

    # Turn into FEniCS Expressions
    variables = [Expression(var, degree=2) for var in variables]

    # Extract variables
    u_e, u_00, u_01, f, g, r, s = variables

    # Define boundary conditions
    boundary_conditions = {0: {'Dirichlet': u_00},   # x = 0
                           1: {'Dirichlet': u_01},   # x = 1
                           2: {'Robin':     (r, s)}, # y = 0
                           3: {'Neumann':   g}}      # y = 1

    # Compute solution
    kappa = Constant(1)
    Nx = Ny = 8
    u = solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
                   degree=1, linear_solver='direct')

    # Compute maximum error at vertices
    mesh = u.function_space().mesh()
    vertex_values_u_e = u_e.compute_vertex_values(mesh)
    vertex_values_u = u.compute_vertex_values(mesh)
    error_max = np.max(np.abs(vertex_values_u_e -
                              vertex_values_u))
    print('error_max =', error_max)

    # Save and plot solution
    vtkfile = File('poisson_extended/solution_bcs.pvd')
    vtkfile << u
    plot(u)

def demo_solvers():
    "Reproduce exact solution to machine precision with different linear solvers"

    # Tolerance for tests
    tol = 1E-10

    # Define exact solution and coefficients
    import sympy as sym
    x, y = sym.symbols('x[0], x[1]')
    u = 1 + x**2 + 2*y**2
    f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)
    f = sym.simplify(f)

    # Generate C/C++ code for UFL expressions
    u_code = sym.printing.ccode(u)
    f_code = sym.printing.ccode(f)

    # Define FEniCS Expressions
    u_e = Expression(u_code, degree=2)
    f = Expression(f_code, degree=2)
    kappa = Constant(1)

    # Define boundary conditions
    boundary_conditions = {0: {'Dirichlet': u_e},
                           1: {'Dirichlet': u_e},
                           2: {'Dirichlet': u_e},
                           3: {'Dirichlet': u_e}}

    # Iterate over meshes and degrees
    for Nx, Ny in [(3, 3), (3, 5), (5, 3), (20, 20)]:
        for degree in 1, 2, 3:
            for linear_solver in 'direct', 'Krylov':
                print('\nSolving on 2 x (%d x %d) mesh with P%d elements '
                      'using solver "%s".' \
                      % (Nx, Ny, degree, linear_solver)),

                # Compute solution
                u = solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
                               degree=degree,
                               linear_solver=linear_solver,
                               abs_tol=0.001*tol,
                               rel_tol=0.001*tol)

                # Compute maximum error at vertices
                mesh = u.function_space().mesh()
                vertex_values_u_e = u_e.compute_vertex_values(mesh)
                vertex_values_u = u.compute_vertex_values(mesh)
                error_max = np.max(np.abs(vertex_values_u_e -
                                          vertex_values_u))
                print('error_max =', error_max)
                assert error_max < tol

if __name__ == '__main__':

    # List of demos
    demos = (demo_test,
             demo_flux,
             demo_convergence_rates,
             demo_structured_mesh,
             demo_bcs,
             demo_solvers)

    # Pick a demo
    for nr in range(len(demos)):
        print('%d: %s (%s)' % (nr, demos[nr].__doc__, demos[nr].__name__))
    print('')
    nr = eval(input('Pick a demo: '))

    # Run demo
    demos[nr]()

    # Hold plot
    import matplotlib.pyplot as plt
    plt.savefig("figure.png")

2 Likes

@dokken: Thanks for your kind reply.

  1. Is the version 2019.1.0, that you mentioned, is the version of dolfin? As per page.3, “Chapter 1: Preliminaries” of FEniCS tutorial, FEniCS uses many components and dolfin being one of them. Please correct me, if my understanding is not right.

  2. How to do I install dolfin 2019.1.0 in my debian machine? I went through a number of websites, like the one below, but I cannot decipher how to install it.
    https://fenicsproject.org/docs/dolfin/dev/python/installation.html#

  3. I am able to download the dolfin 2019.2.0 from these websites:
    https://launchpad.net/debian/+source/dolfin
    https://launchpad.net/debian/+source/dolfin/2019.2.0~git20200218.027d9cc-6

Would you advice how to install dolfin 2019.1.0 or even earlier version, in geneal?
Thanks

  1. To clarify the diference between FEniCS and dolfin. FEniCS is the collection of all components used to solve the finite element problem, this includes UFL, FFC, FIAT and dolfin. Dolfin is the python interface of fenics, and can be imported by either calling from dolfin import * or from fenics import *, as these to commands are equivalent.

  2. You can use the 2019.2.0 version of dolfin (That should be the latest development version of FEniCS).
    In general, I tend to use docker to run dolfin, then you can easily access different versions of dolfin, as images for all version from 2016.1.0 to 2019.2.0 is available.
    A throughout description of installing dolfin/fenics, can be found at: https://fenicsproject.org/download/

  3. I have not used the Launchpad Debian package, but I guess @dparsons might know how to install and use those.

1 Like

If you have 2019.2 available then why use 2019.1? They are code-compatible. Just use the 2019.2 snapshot.

Perhaps you mean how to install 2018.1 rather than 2019.1?

Debian keeps an archive of old builds at http://snapshot.debian.org/ . There may be difficulties with library dependencies, an old version might not install cleanly.

Ubuntu doesn’t have an archive of old versions, they delete superseded builds.

You could try the FEniCS PPA. Click on the “View package details” link on the front page, it takes you to https://launchpad.net/~fenics-packages/+archive/ubuntu/fenics/+packages . The links there provide some builds to older versions for various FEniCS releases (generally, versions newer than the standard Ubuntu build for the given release).

Older versions are not routinely made available though. If you really need them and these links above don’t have the ones you want, you can rebuild the packages yourself.

1 Like

@dparsons: Thanks for response. I never installed an old version by myself. Would you assist me in doing so? Let’s say I would like to install dolfin 2019.2.0
In Debian, /etc/apt/sources.list, I added these 4 lines

# For FEniCS
deb     http://snapshot.debian.org/package/dolfin/2019.2.0~git20200218.027d9cc-6/ buster main
deb-src http://snapshot.debian.org/package/dolfin/2019.2.0~git20200218.027d9cc-6/ buster main
deb     http://snapshot.debian.org/package/dolfin/2019.2.0~git20200218.027d9cc-6/ buster/updates main
deb-src http://snapshot.debian.org/package/dolfin/2019.2.0~git20200218.027d9cc-6/ buster/updates main

Then I used sudo apt-get update, I received the following error messages.

Ign:7 http://snapshot.debian.org/package/dolfin/2019.1.0-10 buster InRelease   
Ign:8 http://snapshot.debian.org/package/dolfin/2019.1.0-10 buster/updates InRelease
Err:9 http://snapshot.debian.org/package/dolfin/2019.1.0-10 buster Release 
  404  Not Found [IP: 2001:630:206:4000:1a1a:0:c13e:ca1b 80]
Err:10 http://snapshot.debian.org/package/dolfin/2019.1.0-10 buster/updates Release
  404  Not Found [IP: 2001:630:206:4000:1a1a:0:c13e:ca1b 80]
Reading package lists... Done
E: The repository 'http://snapshot.debian.org/package/dolfin/2019.1.0-10 buster Release' does not have a Release file.
N: Updating from such a repository can't be done securely, and is therefore disabled by default.
N: See apt-secure(8) manpage for repository creation and user configuration details.
E: The repository 'http://snapshot.debian.org/package/dolfin/2019.1.0-10 buster/updates Release' does not have a Release file.
N: Updating from such a repository can't be done securely, and is therefore disabled by default.
N: See apt-secure(8) manpage for repository creation and user configuration details.

Basically, I do not know how to build a package by myself or even install it. It could be any version. Would you, @dparsons, be kind enough to walk me through? Thanks.

What version of Ubuntu or Debian are you working with? What is written in /etc/debian_version ?

I am using Debian 10.4.

Please answer the rest of the question.

cat /etc/debian_version

I typed cat /etc/debian_version and I received 10.4. My apology is there some other questions you have asked and I have not seen?

Thanks for the clarification. I was expecting it to return a named version like “buster”. 10.4 is ok.

What version of dolfin is provided with Debian 10.4?
(i.e. I am trying to understand why you are trying to install a different version to the one already provided)

1.) I have installed 2019.1 insalled in my machine, most likely using docker.

2.) However, when I typed dolfin-version in terminal outside docker, I received 2018.1.0.

3.) I am not sure which version the Spyder is using, most likely 2018.1.0. I do not know how can I tell Spyder to choose a specific packge to install.

This is to start using the manuals to learn the FEniCS, even if it is older version, not get bogged down and I don’t know if this is related to a recent package update.

As you could see, I tried many ways for the last one month to get started with using FEniCS!

Remove and purge completely all versions installed by docker or pip.
Remove your snapshot entries from /etc/apt/sources.list

Run sudo apt-get update and make sure it is updating against the standard archives for debian 10.4 buster/stable.

Then install FENiCS using

sudo apt-get install fenics

Check the installation is in place with

python3 -c "import dolfin; print(dolfin.__version__)"

1.) I have removed the snapshot entries from the /etc/apt/sources.list.

2.) How to remove and purge completely all versions installed by docker and pip? Frankly I don’t know what is pip. I have seen the word pip before.

Try

pip  list --user

to see if you have modules user-installed using pip. If dolfin or other fenics components are listed there, remove them with

pip uninstall --user dolfin

It sounds like pip is probably not the problem though. You’ve been using docker, not pip.

I received the following after typingpip list --user

Package        Version
-------------- -------
pip-autoremove 0.9.1
WARNING: You are using pip version 20.1; however, version 20.1.1 is available.
You should consider upgrading via the '/usr/bin/python3 -m pip install --upgrade pip' command.

1.) Shall I type the following commands?

apt-get remove fenics
apt-get purge fenics

apt-get install fenics
python3 -c "import dolfin; print(dolfin.__version__)"

2.) I am not sure what you mean by

make sure it is updating against the standard archives for debian 10.4 buster/stable.

3.) I have the following in my sources.list

deb http://deb.debian.org/debian/ buster main non-free contrib
deb-src http://deb.debian.org/debian/ buster main non-free contrib

deb http://security.debian.org/debian-security buster/updates main contrib non-free
deb-src http://security.debian.org/debian-security buster/updates main contrib non-free

# buster-updates, previously known as 'volatile'
deb http://deb.debian.org/debian/ buster-updates main contrib non-free
deb-src http://deb.debian.org/debian/ buster-updates main contrib non-free
deb [arch=amd64] https://download.docker.com/linux/debian buster stable
# deb-src [arch=amd64] https://download.docker.com/linux/debian buster stable

Right, so your sources.list is configured against the buster archive (which is 10.4). You can check it’s being accessed when you run sudo apt-get update.