Mixed space expressions


60
views
0
18 days ago by
I have a mixed space problem that has a nonlinear expression. I'd like the functions that go into the expression to be updated on each iteration of the nonlinear solver such that the functions that go into the expression do not have to be updated on each time so I can use the built-in nonlinear solver. Below I have a rough example of what I'm trying to achieve with a bogus PDE that only serves to demonstrate what I'm trying to do. Is there any variant of the function assigner that will behave as a "copy by reference" (soft copy) assignment rather than a "copy by value"? Specifically in this example, when "prev_sol" gets updated in the nonlinear solver, I'd like for the functions to be updated respectively so that expr gets updated. Thanks!

import sys
import fenics as fem
import numpy as np


def picard_solver(bilin, lin, u, u_1, a, b, c, d):
    eps = 1.0
    tol = 1e-5
    iter = 0
    maxiter = 25
    while eps > tol and iter < maxiter:
        iter += 1
        fem.solve(bilin == lin, u)
        diff = u.vector().get_local() - u_1.vector().get_local()
        eps = np.linalg.norm(diff, ord=np.Inf)
        print('iter={}, norm={}'.format(iter, eps))
        a.assign(u.split(True)[0])
        b.assign(u.split(True)[0])
        c.assign(u.split(True)[0])
        d.assign(u.split(True)[0])
        u_1.assign(u)


def main():
    mesh = fem.IntervalMesh(10, 0, 3)

    P = fem.FiniteElement('CG', mesh.ufl_cell(), degree=1)
    element = fem.MixedElement([P, P, P, P])
    V1 = fem.FunctionSpace(mesh, P)
    V = fem.FunctionSpace(mesh, element)

    prev_sol = fem.Function(V)
    a = fem.Function(V1)
    b = fem.Function(V1)
    c = fem.Function(V1)
    d = fem.Function(V1)

    assigner = fem.FunctionAssigner(V, [V1, V1, V1, V1])
    assigner.assign(prev_sol, [a, b, c, d])

    expr = fem.Expression('a*b*c*d', a=a, b=b, c=c, d=d, degree=1)

    a.vector()[:] = np.linspace(0, 3, 11)[fem.dof_to_vertex_map(V1)]
    b.vector()[:] = a.vector()
    c.vector()[:] = a.vector()
    d.vector()[:] = a.vector()

    # wanted result: [0.  0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3. ]
    # actual result: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
    answer = fem.interpolate(prev_sol.sub(0), V1).vector().get_local()[fem.vertex_to_dof_map(V1)]
    print(answer)

    assigner.assign(prev_sol, [a, b, c, d])

    # result: [0.  0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3. ]
    answer = fem.interpolate(prev_sol.sub(0), V1).vector().get_local()[fem.vertex_to_dof_map(V1)]
    print(answer)

    # random problem:
    u = fem.TrialFunction(V)
    v = fem.TestFunction(V)

    dx = fem.Measure('dx', domain=mesh)
    F = expr*fem.dot(u, v)*dx

    system = fem.Function(V)
    system.assign(prev_sol)
    
    # want to be able to avoid assigning a, b, c, d on each interval
    picard_solver(fem.lhs(F), fem.rhs(F), system, prev_sol, a, b, c, d)

    return 0

if __name__ == '__main__':
    sys.exit(main())​
Community: FEniCS Project
you could try to use an C++ expression where you point directly to the functions you want to be automatically updated... not sure how complicated is that..

Here in the example they show and do something pointing to a CellFunction https://fenicsproject.org/olddocs/dolfin/1.5.0/python/programmers-reference/functions/expression/Expression.html
written 15 days ago by lhdamiani  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »