dolfin team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

On 14 June 2011 23:02, Martin Sandve Alnæs <martinal@xxxxxxxxx> wrote:
> I think one of the main problems with VariationalProblem
> is the mix of problem and solver in one abstraction.
>
> A problem is usually formulated:
>
> Find u in V such that
>  a = L in Omega // pde
>  u = g on Gamma // bcs
>
> Find u in V such that
>  F = 0 in Omega // pde
>  u = g on Gamma // bcs
>
> Here's my suggested variation of VariationalProblem:
>
> # Generic form:
> problem = VariationalProblem(unknown, lhsform, rhsform, bcs)
> solve(problem)
>
> # Examples and variations follow, all consistent with the above:
>
> # Typical linear system
> p = VariationalProblem(u, a, L, bcs)
> # Note that ordering and sign is consistent with a linear system
> p = VariationalProblem(u, J, -F, bcs)
> # Possibly(!) nonlinear system F(u) = 0
> p = VariationalProblem(u, F, 0, bcs)
>
> # Note the 0 which maybe makes this less ambiguous(?)
>
> # The fully automated variant of solve:
> solve(p) # Autodetect linear or nonlinear and differentiate if necessary
>
> # With optional behaviour changes:
> solve(p, nonlinear=True) # Force nonlinear solver
> solve(p, nonlinear=False) # Force linear solver
> solve(p, u=u2) # place solution in u2 instead of p.u, useful for linear solvers
> solve(p, solver=mysolver) # Use nondefault solver (need solver API)
>
> # Pseudocode for solve, showing it's feasible:
> def solve(p, nonlinear=None, solver="auto"):
>    # Allow alternative solvers in some way
>    if solver != "auto":
>        check properties of solver and adjust what is necessary below
>
>    # Differentiate if necessary (skip if we use a solver that doesn't need it):
>    if isinstance(p.rhs, int) and p.rhs == 0:
>        lhs, rhs = derivative(p.lhs, p.u), -p.lhs
>    else:
>        lhs, rhs = p.lhs, p.rhs
>
>    # Detect nonlinear or linear equation unless forced
>    if nonlinear is None: # Set to True or False to force this
>        # I have this function in my private sandbox
>        nonlinear = ufl.algorithms.depends_on(lhs, p.u)
>
>    # Solve with linear or nonlinear default solver
>    if solver == "auto":
>        if nonlinear:
>            return solve_nonlinear(p.u, lhs, rhs, p.bcs)
>        else:
>            return solve_linear(p.u, lhs, rhs, p.bcs)
>    else:
>        return solver.solve(...)
>
>
> Good properties of this solution:
> - Makes the VariationalProblem independent of the solver method
> - Allows explicit or automatic choice of nonlinear/linear solver
> - Consistent ordering of forms given to VariationalProblem
> - Binds u to the problem, which is in accordance with (*) and allows
> differentiation
> There may be something I've forgotten of course, so speak up.
>
> The binding of u to the problem is not really necessary in
> the case of a linear problem, but makes things consistent.
> With solve(p, u=u2) we still allow reuse of the VariationalProblem.
>
> Martin

Plus this is a fairly small iteration from the current design, with
the book and all...

Martin