Efficiently solving Poisson eqn repeatedly

dolfinx
1

Hi all,

For my application, I am solving a Poisson equation repeatedly with a different nonconstant RHS and zero Dirichlet boundary conditions.

If I’m not mistaken, I could simply work out the LU decomposition of the LHS matrix and store it somewhere. Then, I’d apply the decomposition repeatedly to a RHS vector to obtain a solution.

However, I’m not sure how to implement this in dolfinx, below a MWE for what I do now:

from random import random
import dolfinx
import numpy as np
from dolfinx import FunctionSpace, Constant, fem, Function
from dolfinx.fem import create_vector, Form
from mpi4py import MPI
from ufl import TrialFunction, TestFunction, inner, dx
from ufl.operators import grad

N = 2 ** 8
mesh = dolfinx.UnitSquareMesh(MPI.COMM_SELF, N, N, dolfinx.cpp.mesh.CellType.triangle)
V = FunctionSpace(mesh, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
uh = Function(V)


A = inner(grad(u), grad(v)) * dx
F = inner(f, v) * dx
problem = fem.LinearProblem(A, F, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})  
problem.solve()  # Calculate LU of LHS

# Start solving all cases
for index in range(500):
    rand = random()
    f.interpolate(lambda x: np.sin(rand*(x[0] + x[1])))  # Something nontrivial, in reality it's a function of the 
    # previous solutions 
    RHS = inner(f, v) * dx
    form = Form(RHS)
    print("Create")
    vector = create_vector(form)
    problem._b = vector
    problem._L = form
    print("Solve")
    problem.solve()

However, I found this by inspecting the source of dolfinx v0.3.0, and not by following any official resources.

In the KSP documentation (KSP: Linear System Solvers — PETSc 3.14.1 documentation) I can find a mention of my use case, where it mentions successive calls to KSPSolve().

So basically my question is, is this the appropriate method to perform this?

Best whishes,

  • Wouter
2

You can set this explicitly to be sure: KSPSetReusePreconditioner.

1 Like
3

There are several things in your code that is not efficient,
1.

You do not need to recompile the form if f is in the same function space. You can simply create your form outside the loop and update f through interpolation.

This should not be required, as the form should be constant in time, and the coefficients in f will be automatically updated during the solve stage.

2 Likes
4

Thank you for the tips,

I wrote this MWE a bit more simplified than what I’m actually calculating as I thought this would needlessly complicate the example.

In reality I’m solving the weak formulation

\int_{D} \nabla t_i(x)\cdot\nabla v(x) dx = \sum_{j < i}\int_{D} a_{j-i}(x) \nabla t_i(x)\cdot\nabla v(x) dx for some decreasing sequencing of functions a_j(x)

which incorporates the gradient of previously calculated solutions t_i. Would your suggestions still work in this case? I don’t think I can simply re-interpolate my right hand side I’d think.

Thanks in advace,

  • Wouter
5

Thank you for your suggestion, I’ll do this to make sure.

  • Wouter
6

If the variational form is constant as a mathematical expression, for instance L=inner(f(t_i),v)*dx + inner(grad(u(t_(i-1)),grad(v))*dxthis can easily be done with the same form by interpolation into a previous solution functionun` representing u(t_i-1). See for instance:
https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code2.html#variational-form

7

Many Thanks to both @dokken and @nate, It has helped me more than I would have thought!

  • Wouter