Dynamic b vector for linear problem

dolfinx
1

Hello everyone,

I am using dolfinx 0.5.2 and trying to solve a static electromagnetic problem. The problem in question needs dynamic creation of vector b for a linear resolution of type Ax=b.
Currently, a similar code is executed at each iteration:

from dolfinx.fem import (
    Function,
    FunctionSpace,
)
from ufl import (
    TestFunction,
    TrialFunction,
    dot,
    dx,
    grad,
    SpatialCoordinate,
)
from dolfinx.fem.petsc import LinearProblem

V = FunctionSpace(mesh, ("CG", 2))

u = TrialFunctionV)
v = TestFunction(V)

x = SpatialCoordinate(mesh)
a = 1 / (2.0 * np.pi * MU0) * (1 / x[0] * dot(grad(u), grad(v))) * dx
L = J * v * dx

u_h = Function(V)
problem = LinearProblem(a, L, u=u_h, bcs=dirichlet_bcs)
u_h = problem.solve()

Where MU0 is a constant and J is a dolfinx function that varies at each iteration. The idea would be to preliminarily create the matrix A, dynamically create the vector b from J, and then solve the linear system at each iteration in order to speed up calculations.

Thank you in advance for your help.

2

What type is J? You have not specified if it is a constant, a coefficient (dolfinx.fem.Function) or anything else?
For this I would probably not use the LinearProblem.solve, as it re-assembles the matrix at every step.
I would rather use dolfinx.fem.petsc.assemble_matrix and create a PETSc.KSP object that can cache things such as the LU factorization.

3

Thank you for the answer. Yes, J is dolfinx.fem.Function.
Meanwhile, I tried something like

import dolfinx
from dolfinx.fem import (
    Function,
    FunctionSpace,
)
from ufl import (
    TestFunction,
    TrialFunction,
    dot,
    dx,
    grad,
    SpatialCoordinate,
)
import numpy as np

V = FunctionSpace(mesh, ("CG", 2))

u = TrialFunctionV)
v = TestFunction(V)

x = SpatialCoordinate(mesh)
a = dolfinx.fem.forms.form(
      1 / (2.0 * np.pi * MU0) * (1 / x[0] * dot(grad(u), grad(v))) * dx
    )
A = dolfinx.fem.petsc.assemble_matrix(a)
A.assemble()

L = dolfinx.fem.forms.form(J * v * dx)
b = dolfinx.fem.petsc.assemble_vector(L)
b.assemble()

from petsc4py import PETSc
ksp = PETSc.KSP().create()
ksp.setOperators(A)

x = PETSc.Vec()
ksp.solve(b, x)

But I get:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple MacOS to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=59
:
system msg for write_line failure : Bad file descriptor
4

You need to call A.assemble() after using this command.

5

use

x_func = fem.Function(V)
x = x_func.vector

You should also send in the mpi communicator in the ksp.create, see: https://github.com/FEniCS/dolfinx/blob/deac4f1cba124f09302f8bb501777f3e25a827c4/python/test/unit/fem/test_assembler.py#L423
or A known analytical solution — FEniCSx tutorial

6

Thank you so much, seems working!