Equivalence of solve(A, x, b) and solve(a ==L, u)

This is confusing, but the solve interface for forms and la objects are
different. See doc string. To answer you questions use:

   solve(A, x, b, ksp, pc)

for la objects. On could expect an error for the prohibited kwarg
"solver_parameters" for the la object version of solve. Please report a
bug! :)

Johan

On 04/09/2012 11:06 PM, Patrick Farrell wrote:
> New question #193119 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/193119
>
> Hi,
>
> I expected the following two calling sequences to be equivalent (assuming solver_parameters["symmetric"] is False):
>
> solve(a == L, u, bcs, solver_parameters=solver_parameters)
> <=>
> A = assemble(a)
> b = assemble(L)
> [bc.apply(A) for bc in bcs]
> [bc.apply(b) for bc in bcs]
> solve(A, u.vector(), b, solver_parameters=solver_parameters)
>
> However, it appears that's not the case for certain values of "linear_solver" and "preconditioner", which I was surprised by. For example, consider the following script:
>
> from dolfin import *
> import numpy
>
> mesh = Rectangle(0, 0, 2.0, 1.0, 5, 5)
> Q = FunctionSpace(mesh, "DG", 1)
> W = VectorFunctionSpace(mesh, "CG", 2) * FunctionSpace(mesh, "CG", 1)
>
> T_ = Function(Q)
> T_.vector()[:] = numpy.array([1.0, 1.0, 0.5000000000767537, 1.0, 0.910936301182697, 0.5000000000767537, 1.0, 1.0, 0.5000000000000719, 1.0, 0.5000000000767537, 0.5000000000000719, 1.0, 1.0, 0.5, 1.0, 0.5000000000000719, 0.5, 1.0, 1.0, 0.5, 1.0, 0.5, 0.5, 1.0, 0.18169011381620925, 0.16447191930341964, 1.0, 0.5, 0.16447191930341964, 0.910936301182697, 0.5000000000767537, 0.5, 0.910936301182697, 0.8804530826008414, 0.5, 0.5000000000767537, 0.5000000000000719, 0.5, 0.5000000000767537, 0.5, 0.5, 0.5000000000000719, 0.5, 0.5, 0.5000000000000719, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.16447191930341964, 0.14411872829141137, 0.5, 0.5, 0.14411872829141137, 0.8804530826008414, 0.5, 0.5, 0.8804530826008414, 0.8558812717085886, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.14411872829141137, 0.11954691739915857, 0.5, 0.5, 0.11954691739915857, 0.8558812717085886, 0.5, 0.5, 0.8558812717085886, 0.8355280806965801, 0.5, 0.5, 0.5, 0.5, 0
.5, 0.5, 0.5, 0.5, 0.5, 0.49999999999992784, 0.5, 0.5, 0.49999999999992784, 0.5, 0.5, 0.4999999999232465, 0.5, 0.49999999999992784, 0.4999999999232465, 0.5, 0.11954691739915857, 0.08906369881730303, 0.5, 0.4999999999232465, 0.08906369881730303, 0.8355280806965801, 0.5, 0.0, 0.8355280806965801, 0.31830988618379075, 0.0, 0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.49999999999992784, 0.0, 0.5, 0.0, 0.0, 0.49999999999992784, 0.4999999999232465, 0.0, 0.49999999999992784, 0.0, 0.0, 0.4999999999232465, 0.08906369881730303, 0.0, 0.4999999999232465, 0.0, 0.0])
>
> eta = exp(ln(2.0)*(1.0 - triangle.x[1]) )
> f = 1.0e6 * T_ * Constant((0.0, -1.0))
>
> def strain(v):
> return sym(grad(v))
>
> (u, p) = TrialFunctions(W)
> (v, q) = TestFunctions(W)
>
> F = (2.0*eta*inner(strain(u), strain(v))*dx + div(v)*p*dx + div(u)*q*dx + inner(f, v)*dx)
> a = lhs(F)
> L = rhs(F)
>
> bottom = DirichletBC(W.sub(0), (0.0, 0.0), "x[1] == 0.0" )
> top = DirichletBC(W.sub(0).sub(1), 0.0, "x[1] == %g" % 1)
> left = DirichletBC(W.sub(0).sub(0), 0.0, "x[0] == 0.0")
> right = DirichletBC(W.sub(0).sub(0), 0.0, "x[0] == %g" % 2)
> bcs = [bottom, top, left, right]
>
> for (pc, ksp) in [("default", "default"), ("amg", "tfqmr")]:
> u = Function(W)
> solve(a == L, u, bcs, solver_parameters={"preconditioner": pc, "linear_solver": ksp})
>
> u2 = Function(W)
> (A, b) = (assemble(a), assemble(L))
> [bc.apply(A) for bc in bcs]
> [bc.apply(b) for bc in bcs]
> solve(A, u2.vector(), b, solver_parameters={"preconditioner": pc, "linear_solver": ksp})
>
> print "Norm for (%s, %s): " % (pc, ksp)
> print (u2.vector() - u.vector()).norm("l2")
>
> It sets up a linear system, and then attempts the two different approaches, with different values of "preconditioner" and "linear_solver". The output on my system is:
>
> [pef@aislinn:~/src/dolfin-bugs]$ python amg.py
> Solving linear variational problem.
> Norm for (default, default):
> 0.0
> Solving linear variational problem.
> Norm for (amg, tfqmr):
> 4514505.89407
>
> which I wasn't expecting.
>
> If I want to use (amg, tfqmr), how should I use solve(A, x, b) to get the exact same answer as solve(a == L, u)?
>
> Thanks!
>
> Patrick
>
>

Ah! That had me scratching my head for ages. Thanks for the answer :-)

I'll report a bug.