Using the previous step to prime next step

dolfinx
1

Hi, I was wondering whether I can abuse similar solutions (e.g. previous time steps), to speedup the solution of my system.
Below is a minimal example of a pulsating acoustic sphere, where I am probing the response of the domain to different frequencies. Can the solver access the previous solution? How would I have to modify this implementation to get more efficient simulations? Does the solver do that automatically?

Thank you

import dolfinx
import numpy as np
import ufl
import gmsh
from mpi4py import MPI
from dolfinx.fem.petsc import LinearProblem

# mesh
gmsh.initialize()
gmsh.model.add("Pulsating-Sphere")

box = gmsh.model.occ.add_rectangle(0.0, 0.0, 0.0, 1., 1.)
circ = gmsh.model.occ.add_disk(0.0, 0.0, 0.0, 0.2, 0.2, )
gmsh.model.occ.cut([(2, box)], [(2, circ)])

gmsh.model.occ.synchronize()
surfs = gmsh.model.occ.get_entities(2)
gmsh.model.add_physical_group(2, [s[1] for s in surfs])
lines = gmsh.model.occ.get_entities(1)
for line in lines:
    com = np.array(gmsh.model.occ.get_center_of_mass(*line))
    if not np.greater(com[:2], np.zeros(2, )).all():
        continue
    if not np.less(com[:2], np.ones(2, )).all():
        continue
    gmsh.model.add_physical_group(1, [line[1]], tag=42)
gmsh.model.add_physical_group(1, [l[1] for l in lines])
gmsh.option.setNumber("Mesh.MeshSizeMax", 343. / 500. / 15.)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)

msh, ct, ft = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, MPI.COMM_SELF, 0
)
ds = ufl.Measure("ds", msh, subdomain_data=ft)

# variational problem
v = dolfinx.fem.FunctionSpace(
    msh, ufl.FiniteElement("Lagrange", msh.ufl_cell(), 2)
)
p = ufl.TrialFunction(v)
xi = ufl.TestFunction(v)
p_sol = dolfinx.fem.Function(v)
p_sol.name = "p"

# physics constants
v0 = 1e-3
s = 1j * 1.2 * 343.

# solve different systems
for f in np.linspace(450., 500., 50):
    k = f * 2. * np.pi / 343.
    ks = k ** 2

    lhs = (
            ufl.inner(ufl.grad(p), ufl.grad(xi)) * ufl.dx
            - (ks * ufl.inner(p, xi) * ufl.dx)
            - (s * k * 0.5 * (1 + 1j) * ufl.inner(p, xi) * ufl.ds)
    )
    rhs = s * k * ufl.inner(v0, xi) * ds(42)

    problem = LinearProblem(
        lhs,
        rhs,
        u=p_sol,
        petsc_options={
            "ksp_type": "preonly",
            "pc_type": "cholesky",
            "pc_factor_mat_solver_type": "mumps",
        },
    )
    problem.solve()
2

At the first iteration in the loop over f, the linear problem is solved, and the solution is stored in p_sol. At the second iteration, the solution currently stored in p_sol is not deleted, and it is passed to the linear solver [*]. After the linear solve is completed, the content of p_sol is overwritten.

[*] then it’s up to to solver to decide what to do with that. In your case, since mumps is a sparse direct solver, I doubt that the current content of p_sol is used at all.

1 Like