Use Previous Solution as Initial Guess

General
, ,
1

Hi. I have a code written using Legacy Dolfin to solve the Green’s function of the Stokes flow. Actually I solve for a force dipole, defining two point forces instead of one, as you see in the code. What I want to do is to use the solution for “u” and “p” from my code and incorporate it as initial guess for my next iteration of the same problem so the solution converges quicker.
Let me explain myself better. In my simulations the center of mass of the dipole shown in the uploaded image moves tiny increments following certain equation of motion. What I want to do is to compute the velocity field generated by the dipole using the code uploaded before, and use this solution as initial guess at the next position calculation. In theory that should speed up things.
I don’t know if I have to migrate my code to fenicsx to do this. I’m new to the Fenics project, so if I do have to migrate to fenicsx could you please help me with this too. Thanks. Here is my code:

from dolfin import *

mesh = UnitCubeMesh(10, 10, 10)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, MixedElement([V.ufl_element(), Q.ufl_element()]))

# Boundaries
def all_boundaries(x, on_boundary):
return x[0] > (1.0 - DOLFIN_EPS) or x[0] < (DOLFIN_EPS) or x[1] > (1.0 - DOLFIN_EPS) or x[1] < (DOLFIN_EPS) or x[2] > (1.0 - DOLFIN_EPS) or x[2] < (DOLFIN_EPS)
def top_y(x, on_boundary): return x[1] > (1.0 - DOLFIN_EPS)

noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, all_boundaries)
bcs = [bc0]

# lid = Constant((1.0, 0.0, 0.0))
# bc1 = DirichletBC(W.sub(0), lid, top_y)
# bcs = [bc0, bc1]

# Define the directional Point Sources
point0 = Point(0.45, 0.5, 0.5)
delta0 = PointSource(W.sub(0).sub(0), point0, -1.0)

point1 = Point(0.55, 0.5, 0.5)
delta1 = PointSource(W.sub(0).sub(0), point1, 1.0)

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))

a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

b = inner(grad(u), grad(v))*dx + p*q*dx

A, bb = assemble_system(a, L, bcs)

# Apply point sources
delta0.apply(bb)
delta1.apply(bb)

P, btmp = assemble_system(b, L, bcs)

solver = KrylovSolver("minres", "amg")
solver.set_operators(A, P)

U = Function(W)
solver.solve(U.vector(), bb)

u, p = U.split()
2

At the end of your code you have U. Just define the residual of the nonlinear problem using that solution U. Just note that in defining the residual you must be using u, p = split(U), instead of u, p = U.split() which you are currently using (e.g., to make a plot). Search the forum using split as keyword for further topics explaining the difference.

1 Like