Hi, I have a problem where I need to solve a similar PDE a lot of times. Consider the following “minimal” example:
Let us say that we have solved the PDE problem of the general FEniCSx tutorial \begin{align} -\nabla^2 u(\mathbf{x}) &= f(\mathbf{x})&&\mathbf{x} \in \Omega\\ u(\mathbf{x}) &= u_D(\mathbf{x})&& \mathbf{x} \in \partial\Omega\\ f(\mathbf{x}) & = -6\\ u_D(\mathbf{x})&=1+x^2+2y^2 \\ \Omega&=[0,1]\times [0,1] \end{align}
Say now that I want to repeat the above calculation in the same geometry and with the same finite elements but with different source and different boundary conditions (still of the Dirichlet type though).
Now, since, if I am not terribly wrong, only the load vector L and the bc object change and not the inner product a is it possible to find the solution of the new problem without a full calculation, i.e. by the action of the inverse of the stiffness matrix that was calculate when the initial PDE was solved?
Also, is the same possible for more complex problems like the static Navier Stokes, for example?
The minimal code follows.
from mpi4py import MPI
import dolfinx
import matplotlib.pyplot as plt
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.FunctionSpace(domain, ("CG", 1))
uD = dolfinx.fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, domain.topology.dim-1, boundary_facets)
dof_coordinates = V.tabulate_dof_coordinates()[boundary_dofs]
bc = dolfinx.fem.dirichletbc(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = dolfinx.fem.Constant(domain, ScalarType(-6))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])
uh = problem.solve()
I apologise for asking something that seems naive.