Computing gradient on a mixed subspace

dolfin-adjoint
1

Hello all,

If I create a function from a mixed space and I want to define two new functions from the linear combinations of the subfunctions, how do I do this? I want to be able to compute the sensitivity with respect to the original function. I tried:

element = FiniteElement("CG", mesh.ufl_cell(), 1) 
mixed_element = MixedElement([element,element]) 
V = FunctionSpace(mesh,mixed_element)
func = Function(V)

singleVspace = V.sub(0).collapse()
Wp = Function(singleVspace)
Wm = Function(singleVspace)
assigner0 = FunctionAssigner(singleVspace, V.sub(0))
assigner1 = FunctionAssigner(singleVspace, V.sub(1))

assigner0.assign(Wp, func.sub(0))
assigner1.assign(Wm, func.sub(1))

Wa = Wp - Wm
Wb = Wp + Wm

but after trying to compute a gradient with respect to the function func, i get an error:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to create function.
*** Reason:  Cannot be created from subspace. Consider collapsing the function space.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  8060a84e8ada0df018772497534092727b06bd56
*** -------------------------------------------------------------------------

Any help is appreciated!

2

Please add a minimal reproducible example, as it is very hard to see where you error is when one cannot reproduce it.

3

Sorry about that, here is a MWE:

from dolfin import *
from dolfin_adjoint import *
import pyadjoint
from ufl import ln, diff
import numpy as np
import os

# Parameters
N_steps = 50
Rg = 1.0
f = 0.4
χN = 14.4


mesh = UnitSquareMesh(32, 32)
print(mesh.geometry().dim())

element = FiniteElement("CG", mesh.ufl_cell(), 1) 
mixed_element = MixedElement([element,element]) 

V = FunctionSpace(mesh,mixed_element) #Mixed element function space
dx = Measure("dx", domain=mesh)

# Initialize fields
InteractionFields = Function(V)
InteractionFields.vector()[:] = np.random.uniform(low=-1.0, high=1.0, size=InteractionFields.vector().size())

Volume = assemble(1.0*dx)

def diff_bilinear_form(q, ϕ, R, W):
    return (-inner(R**2*grad(ϕ), grad(q)) - ϕ*W*q)*dx

def Crank_Nicholson_residual(q, q_prev, ϕ, R, W, ΔS):
    return inner(q, ϕ)*dx - inner(q_prev, ϕ)*dx - 0.5*ΔS*(diff_bilinear_form(q, ϕ, R, W) + diff_bilinear_form(q_prev, ϕ, R, W))

def assemble_forms(q, ϕ, q_prev, R, ΔS, W):
    Q_form = Crank_Nicholson_residual(q, q_prev, ϕ, R, W, ΔS)
    aq, Lq = lhs(Q_form), rhs(Q_form)
    Aq = assemble(aq)
    
    solverQ = PETScKrylovSolver("cg", "hypre_amg")
    solverQ.set_operator(Aq)

    return solverQ, Lq

def UpdateQ(s,f,solverA,solverB,La,Lb,sol):
    if s < f:
            b = assemble(La)
            solverA.solve(sol.vector(), b)
    else:
            b = assemble(Lb)
            solverB.solve(sol.vector(), b)

def PartitionFunction(N_steps, R, InteractionFields, f):
    ΔS = 1.0/N_steps
    
    # Density functions

    singleVspace = V.sub(0).collapse()
    # Solutions
    qsol = Function(singleVspace)

    Wp = Function(singleVspace)
    Wm = Function(singleVspace)

    assigner0 = FunctionAssigner(singleVspace, V.sub(0))
    assigner1 = FunctionAssigner(singleVspace, V.sub(1))

    assigner0.assign(Wp,InteractionFields.sub(0))
    assigner1.assign(Wm,InteractionFields.sub(1))

    Wa = Wp - Wm
    Wb = Wp + Wm

    # Set initial condition
    qsol.vector()[:] = 1.0

    q = TrialFunction(singleVspace)
    ϕ = TestFunction(singleVspace)

    solverAq, Laq = assemble_forms(q,ϕ,qsol,R,ΔS,Wa)
    solverBq, Lbq = assemble_forms(q,ϕ,qsol,R,ΔS,Wb)

    for i in range(N_steps):
        s = (i+1)*ΔS
        UpdateQ(s,f,solverAq,solverBq,Laq,Lbq,qsol) #Update propagator q

    Qq = assemble(qsol*dx)

    return Qq


Qq  = PartitionFunction(N_steps, Rg, InteractionFields, f)
controlFields = Control(InteractionFields)
δQδW = compute_gradient(Qq, controlFields)
1 Like
4

You can rewrite your code in the following way:

from dolfin import *
from dolfin_adjoint import *
import numpy as np

# Parameters
N_steps = 50
Rg = 1.0
f = 0.4
χN = 14.4


mesh = UnitSquareMesh(32, 32)
print(mesh.geometry().dim())

element = FiniteElement("CG", mesh.ufl_cell(), 1)
mixed_element = MixedElement([element, element])

V = FunctionSpace(mesh, mixed_element)  # Mixed element function space
dx = Measure("dx", domain=mesh)

# Initialize fields
InteractionFields = Function(V)
InteractionFields.vector()[:] = np.random.uniform(
    low=-1.0, high=1.0, size=InteractionFields.vector().size()
)

Volume = assemble(1.0 * dx)


def diff_bilinear_form(q, ϕ, R, W):
    return (-inner(R**2 * grad(ϕ), grad(q)) - ϕ * W * q) * dx


def Crank_Nicholson_residual(q, q_prev, ϕ, R, W, ΔS):
    return (
        inner(q, ϕ) * dx
        - inner(q_prev, ϕ) * dx
        - 0.5
        * ΔS
        * (diff_bilinear_form(q, ϕ, R, W) + diff_bilinear_form(q_prev, ϕ, R, W))
    )


def assemble_forms(q, ϕ, q_prev, R, ΔS, W):
    Q_form = Crank_Nicholson_residual(q, q_prev, ϕ, R, W, ΔS)
    aq, Lq = lhs(Q_form), rhs(Q_form)
    Aq = assemble(aq)

    solverQ = PETScKrylovSolver("cg", "hypre_amg")
    solverQ.set_operator(Aq)

    return solverQ, Lq


def UpdateQ(s, f, solverA, solverB, La, Lb, sol):
    if s < f:
        b = assemble(La)
        solverA.solve(sol.vector(), b)
    else:
        b = assemble(Lb)
        solverB.solve(sol.vector(), b)


def PartitionFunction(N_steps, R, InteractionFields, f):
    ΔS = 1.0 / N_steps

    # Density functions

    V0 = V.sub(0).collapse()
    V1 = V.sub(1).collapse()
    # Solutions
    qsol = Function(V0)

    Wp = Function(V0)
    Wm = Function(V1)

    assigner = FunctionAssigner([V0, V1], V)

    assigner.assign([Wp, Wm], InteractionFields)

    Wa = Wp - Wm
    Wb = Wp + Wm

    # Set initial condition
    qsol.vector()[:] = 1.0

    q = TrialFunction(V0)
    ϕ = TestFunction(V0)

    solverAq, Laq = assemble_forms(q, ϕ, qsol, R, ΔS, Wa)
    solverBq, Lbq = assemble_forms(q, ϕ, qsol, R, ΔS, Wb)

    for i in range(N_steps):
        s = (i + 1) * ΔS
        UpdateQ(s, f, solverAq, solverBq, Laq, Lbq, qsol)  # Update propagator q

    Qq = assemble(qsol * dx)

    return Qq


Qq = PartitionFunction(N_steps, Rg, InteractionFields, f)
controlFields = Control(InteractionFields)
δQδW = compute_gradient(Qq, controlFields)
5

Thanks a lot , this now runs smoothly! Is it possible to extract the gradients with respect to the individual components from the gradient of the mixed element?

6

I am not sure what you are asking for here? Would you like to only get the gradient with respect to a sub-component of

?
Then I would start by trying to set Control(InteractionFields.sub(0)) and see if that works.