Mutliphysics Iterative Solutions

[B-Nazari]

Hi,
Is it possible to solve a coupled problem in multistage iterative manner? Let's say if the forcing for problem 1 is a function of solution of problem 2, can solve calls be nested so that after solving problem 1 the solution is used to compute forcing for problem 2 which will be used for next time step of problem 1 and so on?
If there's already an example usage I would appreciate if someone shares the link. Also I know adjoint modeling is possible but that requires changing mathematical formulation and may even be more time consuming in the inverse stage where a loose coupling might be enough.
Thanks

[colinjcotter]

Anything that you can write as a series of solves of problems that can be written in weak form, can be done in Firedrake.

[B-Nazari]

Thanks, I will try to look into it to see if I can rephrase all the equations in weak form.

[rckirby]

Your proposed algorithm might be possible to express via PETSc options without writing separate problems. Write your problem as a monolithic, fully-coupled system for a given time-step. Then, the fractional method step you describe is probably going to arise by letting the nonlinear solve be ksponly and the linear solve as preonly. The preconditioner is then a multiplicative field split (solve the first block of equations for the first unknown, stuff the result into the second block), and you invert the diagonal blocks however you like.

This has the advantage of letting you switch between monolithic and lagged solvers to increase stability/accuracy as desired without rewriting an application code. Adjoints should work in either case as long as the solver options are respected under adjoints.

[B-Nazari]

Thanks. This sounds interesting I need to read more about PETSC options.

[B-Nazari]

So it would be something like this?

`from firedrake import *

mesh = UnitSquareMesh(10, 10)

V = FunctionSpace(mesh, 'P', 1)
W = FunctionSpace(mesh, 'P', 1)

u, v = TrialFunction(V), TrialFunction(W)
phi, psi = TestFunction(V), TestFunction(W)

u_sol = Function(V)
v_sol = Function(W)

bc1 = DirichletBC(V, 0, 'on_boundary')
bc2 = DirichletBC(W, 0, 'on_boundary')

def g(v):
return 10 * sin(v)

a1 = dot(grad(u), grad(phi)) * dx
L1 = g(v_sol) * phi * dx

a2 = v * psi * dx + dot(grad(v), grad(psi)) * dx
L2 = Constant(1) * psi * dx

a = a1 + a2
L = L1 + L2

solve(a == L, [u_sol, v_sol], bcs=[bc1, bc2], solver_parameters={
"snes_type": "ksponly",
"ksp_type": "preonly",
"pc_type": "fieldsplit",
"pc_fieldsplit_type": "multiplicative",
})`