Time-dependent Navier Stokes flow control

dolfin-adjoint
1

I’m trying to control the recirculation area after a backward facing step using dolfin_adjoint.
The control is performed wrt the amplitude of a inflow parabolic profile. I’ve implemented the forward solver as an incremental Chorin Teman method. The simulation crashes when

    dJd0 = compute_gradient(J, control)

is called and the error message displays

Traceback (most recent call last):
  File "non_time_control.py", line 393, in <module>
    optimize()
  File "non_time_control.py", line 361, in optimize
    dJd0 = compute_gradient(J, control)
  File "/home/fenics/.local/lib/python3.6/site-packages/pyadjoint/drivers.py", line 29, in compute_gradient
    tape.evaluate_adj(markings=True)
  File "/home/fenics/.local/lib/python3.6/site-packages/pyadjoint/tape.py", line 140, in evaluate_adj
    self._blocks[i].evaluate_adj(markings=markings)
  File "/home/fenics/.local/lib/python3.6/site-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
  File "/home/fenics/.local/lib/python3.6/site-packages/pyadjoint/block.py", line 133, in evaluate_adj
    prepared)
  File "/home/fenics/.local/lib/python3.6/site-packages/fenics_adjoint/types/dirichletbc.py", line 148, in evaluate_adj_component
    adj_output = r
UnboundLocalError: local variable 'r' referenced before assignment

I am not able to understand how to solve it. Do you have any suggestion?
Thank you

from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *
import numpy as np
import os, sys
set_log_level(LogLevel.ERROR)


#Class for the bd term
class Source(UserExpression):
    def __init__(self, omega=Constant(1e-1), **kwargs):
        """ Construct the source function """
        super().__init__(self,**kwargs)
        self.t = 0.0
        self.omega = omega

    def eval(self, value, x):
        """ Evaluate the source function """
        
        value[1] = float(self.omega)*(x[0]-0.95)*(1-x[0])/0.0025
        value[0] = Constant((0.0))
    

    def value_shape(self):
        return (2,)


class SourceDerivative(UserExpression):
    def __init__(self, omega=Constant(1e-1), Source=None, **kwargs):
        """ Construct the source function derivative """
        super().__init__(**kwargs)
        self.t = 0.0
        self.omega = omega
        self.source = Source  # needed to get the matching time instant

    def eval(self, value, x):
        """ Evaluate the source function's derivative """
        
        value[1] = (x[0]-0.95)*(1-x[0])/0.0025
        value[0] = Constant((0.0))
    
    def value_shape(self):
        return (2,)



# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p, mu):
    return 2*mu*epsilon(u) - p*Identity(len(u))


def inflow_profile(mesh,dg):
    bot = mesh.coordinates().min(axis=0)[1]+0.1
    top = mesh.coordinates().max(axis=0)[1]

    # width of inlet channel
    H = top - bot

    # # 
    U_in = 1.5
    return Expression(('-4*U_in*(x[1]-bot)*(x[1]-top)/H/H' , '0'), bot=bot, top=top, H=H, U_in=U_in, degree=dg)



def forward(g, annotate=False):

    #Importing mesh
    mesh = Mesh()
    f    = HDF5File(mesh.mpi_comm(), 'meshes/our_mesh.h5','r')
    f.read(mesh, 'mesh', False)

    # function that extracts the entities of dimension equal to the one of the 
    # topology minus one from a mesh
    surfaces = MeshFunction('size_t', mesh, mesh.topology().dim()-1)     

    # apply mesh function to extract facets from mesh and store them into surfaces
    f.read(surfaces, 'facet')                                           


    # These boundary condition tags should be hardcoded by gmsh during generation
    inlet_tag = 1
    outlet_tag = 2
    wall_tag1 = 3 
    # bottom before jet
    wall_tag2 = 4
    # step and bottom after step
    wall_tag3 = 5 
    # top wall
    control_tag = 6           
    
    V_h = VectorElement("CG", mesh.ufl_cell(), 2)
    Q_h = FiniteElement("CG", mesh.ufl_cell(), 1)
    W = FunctionSpace(mesh, V_h * Q_h)
    V_init, Q_init = W.split()

    v, q = TestFunctions(W)
    x = TrialFunction(W)
    u, p = split(x)
    s = Function(W, name="State")
    V_collapse = V_init.collapse()
    
    
    nu = 1E-3

    # Our functional requires the computation of a boundary integral
    # over :math:`\partial \Omega_{\textrm{circle}}`.  Therefore, we need
    # to create a measure for this integral, which will be accessible as
    # :py:data:`ds(2)` in the definition of the functional. In addition, we
    # define our strong Dirichlet boundary conditions.



    # # Define boundary conditions
    # # Inlet velocity
    bcu_inlet = DirichletBC(W.sub(0), inflow_profile(mesh,2),surfaces,inlet_tag) 
    # # No slip
    bcu_wall1 = DirichletBC(W.sub(0), Constant((0., 0.)),surfaces, wall_tag1)
    bcu_wall2 = DirichletBC(W.sub(0), Constant((0., 0.)),surfaces, wall_tag2)
    bcu_wall3 = DirichletBC(W.sub(0), Constant((0., 0.)),surfaces, wall_tag3)
        
    # # Fixing outflow pressure
    bcp_outflow = DirichletBC(W.sub(1), Constant(0.),surfaces, outlet_tag)  

    bcu_jet = DirichletBC(W.sub(0), g, surfaces ,  control_tag)

    bcu = [bcu_inlet,bcu_wall1,bcu_wall2,bcu_wall3,bcu_jet]

    bcp = [bcp_outflow]

    bcs = [bcu_inlet,bcu_wall1,bcu_wall2,bcu_wall3,bcu_jet, bcp_outflow]


    a = (nu*inner(grad(u), grad(v))*dx
        - inner(p, div(v))*dx
        - inner(q, div(u))*dx
        )
    L = inner(Constant((0, 0)), v)*dx

    # Next we assemble and solve the system once to record it with
    # :py:mod:`dolin-adjoint`.

    A, b = assemble_system(a, L, bcs)
    solve(A, s.vector(), b)

    # Next we define the functional of interest :math:`J`, the
    # optimisation parameter :math:`g`, and create the reduced
    # functional.

    u_init, p_init = split(s)
    
    
    
    #Final time and time dicretization
    T = 0.5
    dtn = 0.0005
    num_steps = int(T/dtn)

    #Flow constants
    mu = 1E-3
    rho = 1  # density
    Re=70 #Reynold's number



    #Defining the functional spaces
    Q = FunctionSpace(mesh, 'CG', 1)
    V = VectorFunctionSpace(mesh, 'CG', 2)

    
    # Define boundary conditions
    # Inlet velocity
    bcu_inlet = DirichletBC(V, inflow_profile(mesh,2),surfaces,inlet_tag) 
    # No slip
    bcu_wall1 = DirichletBC(V, Constant((0., 0.)),surfaces, wall_tag1)
    bcu_wall2 = DirichletBC(V, Constant((0., 0.)),surfaces, wall_tag2)
    bcu_wall3 = DirichletBC(V, Constant((0., 0.)),surfaces, wall_tag3)
        
    # Fixing outflow pressure
    bcp_outflow = DirichletBC(Q, Constant(0.),surfaces, outlet_tag)  

    bcu_jet = DirichletBC(V, g, surfaces, control_tag)

    bcu = [bcu_inlet,bcu_wall1,bcu_wall2,bcu_wall3,bcu_jet]

    bcp = [bcp_outflow]

    #Define trial and test functions
    u,v = TrialFunction(V), TestFunction(V)
    p,q =TrialFunction(Q), TestFunction(Q)
    
    # Define functions for solutions at previous and current time steps
    u_n = Function(V,annotate=annotate)
    u_ = Function(V,annotate=annotate)
    p_n = Function(Q,annotate=annotate)
    p_ = Function(Q,annotate=annotate)
    u_aux = Function(V)

    
    # Define expressions used in variational forms
    U = 0.5*(u_n + u)
    n = FacetNormal(mesh)
    f = Constant((0, 0))
    ix= Constant((1,0))
    k = Constant(dtn)
    mu = Constant(mu)

    # Define variational problem for step 1

    F1 = rho*dot((u - u_n) / k, v)*dx \
        + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
        + inner(sigma(U, p_n,mu), epsilon(v))*dx \
        + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
        - dot(f, v)*dx
    a1 = lhs(F1)
    L1 = rhs(F1)
    # Define variational problem for step 2
    a2 = dot(nabla_grad(p), nabla_grad(q))*dx
    L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx
    # Define variational problem for step 3
    a3 = dot(u, v)*dx
    L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
    # Assemble matrices
    A1 = assemble(a1)
    A2 = assemble(a2)
    A3 = assemble(a3)
    # Apply boundary conditions to matrices
    [bc.apply(A1) for bc in bcu]
    [bc.apply(A2) for bc in bcp]
    [bc.apply(A3) for bc in bcu]

    As = [A1,A2,A3]

    # Create solvers
    solvers = list(map(lambda x: LUSolver(), range(3)))


    # Set matrices
    for s, A in zip(solvers, As):
        s.set_operator(A)




    u_start = project(u_init,V)
    p_start = project(p_init,Q)

    u_.assign(u_start,annotate=annotate)
    p_.assign(p_start,annotate=annotate)




    i = 1
    t = 0.0  # Initial time
    times = [t, ]


    import tqdm
    timer = tqdm.tqdm("Solving forward problem", total=num_steps)
    outfile_u = XDMFFile("output/u.xdmf")
    outfile_p = XDMFFile("output/p.xdmf")

    outfile_u.write(u_,t)
    outfile_p.write(p_,t)

    Jtemp = assemble(inner(curl(u_),curl(u_))*dx)
    Jlist = [Jtemp]

        
    while t <= T:
        
        timer.update(1)
        # Step 1: Tentative velocity step
        b1 = assemble(L1)
        [bc.apply(b1) for bc in bcu]
        solvers[0].solve(u_.vector(), b1,annotate=annotate)
        # Step 2: Pressure correction step
        b2 = assemble(L2)
        [bc.apply(b2) for bc in bcp]
        solvers[1].solve(p_.vector(), b2,annotate=annotate)
        # Step 3: Velocity correction step
        b3 = assemble(L3)
        solvers[2].solve(u_.vector(), b3,annotate=annotate)
        
        # Plot solution
        # print("Current value of u_ at time", t, "is:", u_.vector().get_local())
        
        # Update previous solution
        u_n.assign(u_,annotate=annotate)
        p_n.assign(p_,annotate=annotate)
        outfile_u.write(u_,t)
        outfile_p.write(p_,t)


        Jtemp = assemble(inner(curl(u_),curl(u_))*dx)
        Jlist.append(Jtemp)
    
        # # Update time
        t += float(dtn)
        

    outfile_u.close()
    outfile_p.close()

    
    
    return u_,Jlist,dtn


# Callback function for the optimizer
# Writes intermediate results to a logfile
def eval_cb(j, m):
    """ The callback function keeping a log """

    print("omega = %15.10e " % float(m))
    print("Functional = %15.10e " % j)


def optimize(dbg=False):
    """ The optimization routine """

    #Final time and time dicretization
    T = 0.5
    dtn = 0.0005
    num_steps = int(T/dtn)

    # Define the control
    Omega = Constant(2e-1)
    source = Source(Omega, degree=2, name="source")

    # provide the coefficient on which this expression depends and its derivative
    source.dependencies = [Omega]
    source.user_defined_derivatives = {Omega: SourceDerivative(Omega, Source=source, degree=2)}
    


    # Execute first time to annotate the tape
    u,Jlist,dtn = forward(source, True)


    # Define the control
    control = Control(Omega)

    J = 0
    for i in range(1, len(Jlist)):
        J += 0.5*(Jlist[i-1] + Jlist[i])*float(dtn)

    # compute the gradient
    dJd0 = compute_gradient(J, control)
    print("gradient = ", float(dJd0))

    h = Constant(0.1)



    # Prepare the reduced functional
    reduced_functional = ReducedFunctional(J, control, eval_cb_post=eval_cb)
    
    #Taylor test for convergence rate
    conv_rate = taylor_test(reduced_functional, Omega, h)

    # Run the optimisation
    omega_opt = minimize(reduced_functional, method="L-BFGS-B",
                         tol=1.0e-12, options={"disp": True, "gtol": 1.0e-12})

    # Print the obtained optimal value for the controls
    print("omega = %f" % float(omega_opt))
 
    return omega_opt

optimize()
2

As you havent provided the mesh, I cannot reproduce your error.

Please also state what version of FEniCS and dolfin-adjoint that you are running, and how you installed them.

3

Fenics version is: 2019.1.0, taken from a docker image and then I installed inside the container dolfin-adjoint using “pip install”.
You can find the mesh at: our_mesh.h5 - Google Drive
Thank you for your time

4

Please note that your problem has alot of time steps. I would suggest you try to reduce the example to say two time steps for now to make it easier for people to debug.

I would simply implement your problem without having custom derivatives:

from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *
import numpy as np
import os
import sys
set_log_level(LogLevel.ERROR)


# Class for the bd term
class Source(UserExpression):
    def __init__(self, omega=Constant(1e-1), **kwargs):
        """ Construct the source function """
        super().__init__(self, **kwargs)
        self.t = 0.0

    def eval(self, value, x):
        """ Evaluate the source function """

        value[1] = (x[0]-0.95)*(1-x[0])/0.0025
        value[0] = Constant((0.0))

    def value_shape(self):
        return (2,)


class SourceDerivative(UserExpression):
    def __init__(self, Source=None, **kwargs):
        """ Construct the source function derivative """
        super().__init__(**kwargs)
        self.t = 0.0
        self.source = Source  # needed to get the matching time instant

    def eval(self, value, x):
        """ Evaluate the source function's derivative """

        value[1] = (x[0]-0.95)*(1-x[0])/0.0025
        value[0] = Constant((0.0))

    def value_shape(self):
        return (2,)


# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor


def sigma(u, p, mu):
    return 2*mu*epsilon(u) - p*Identity(len(u))


def inflow_profile(mesh, dg):
    bot = mesh.coordinates().min(axis=0)[1]+0.1
    top = mesh.coordinates().max(axis=0)[1]

    # width of inlet channel
    H = top - bot

    # #
    U_in = 1.5
    return Expression(('-4*U_in*(x[1]-bot)*(x[1]-top)/H/H', '0'), bot=bot, top=top, H=H, U_in=U_in, degree=dg)


def forward(g, omega, annotate=False):

    # Importing mesh
    mesh = Mesh()
    f = HDF5File(mesh.mpi_comm(), 'our_mesh.h5', 'r')
    f.read(mesh, 'mesh', False)

    # function that extracts the entities of dimension equal to the one of the
    # topology minus one from a mesh
    surfaces = MeshFunction('size_t', mesh, mesh.topology().dim()-1)

    # apply mesh function to extract facets from mesh and store them into surfaces
    f.read(surfaces, 'facet')

    # These boundary condition tags should be hardcoded by gmsh during generation
    inlet_tag = 1
    outlet_tag = 2
    wall_tag1 = 3
    # bottom before jet
    wall_tag2 = 4
    # step and bottom after step
    wall_tag3 = 5
    # top wall
    control_tag = 6

    V_h = VectorElement("CG", mesh.ufl_cell(), 2)
    Q_h = FiniteElement("CG", mesh.ufl_cell(), 1)
    W = FunctionSpace(mesh, V_h * Q_h)
    V_init, Q_init = W.split()

    v, q = TestFunctions(W)
    x = TrialFunction(W)
    u, p = split(x)
    s = Function(W, name="State")
    V_collapse = V_init.collapse()

    nu = 1E-3

    # Our functional requires the computation of a boundary integral
    # over :math:`\partial \Omega_{\textrm{circle}}`.  Therefore, we need
    # to create a measure for this integral, which will be accessible as
    # :py:data:`ds(2)` in the definition of the functional. In addition, we
    # define our strong Dirichlet boundary conditions.

    # # Define boundary conditions
    # # Inlet velocity
    bcu_inlet = DirichletBC(
        W.sub(0), inflow_profile(mesh, 2), surfaces, inlet_tag)
    # # No slip
    bcu_wall1 = DirichletBC(W.sub(0), Constant((0., 0.)), surfaces, wall_tag1)
    bcu_wall2 = DirichletBC(W.sub(0), Constant((0., 0.)), surfaces, wall_tag2)
    bcu_wall3 = DirichletBC(W.sub(0), Constant((0., 0.)), surfaces, wall_tag3)

    # # Fixing outflow pressure
    bcp_outflow = DirichletBC(W.sub(1), Constant(0.), surfaces, outlet_tag)

    bcu_jet = DirichletBC(W.sub(0), g, surfaces,  control_tag)

    bcu = [bcu_inlet, bcu_wall1, bcu_wall2, bcu_wall3, bcu_jet]

    bcp = [bcp_outflow]

    bcs = [bcu_inlet, bcu_wall1, bcu_wall2, bcu_wall3, bcu_jet, bcp_outflow]

    a = (nu*inner(grad(u), grad(v))*dx
         - inner(p, div(v))*dx
         - inner(q, div(u))*dx
         )
    L = inner(Constant((0, 0)), v)*dx

    # Next we assemble and solve the system once to record it with
    # :py:mod:`dolin-adjoint`.

    A, b = assemble_system(a, L, bcs)
    solve(A, s.vector(), b)

    # Next we define the functional of interest :math:`J`, the
    # optimisation parameter :math:`g`, and create the reduced
    # functional.

    u_init, p_init = split(s)

    # Final time and time dicretization
    T = 0.5
    dtn = 0.1
    num_steps = int(T/dtn)

    # Flow constants
    mu = 1E-3
    rho = 1  # density
    Re = 70  # Reynold's number

    # Defining the functional spaces
    Q = FunctionSpace(mesh, 'CG', 1)
    V = VectorFunctionSpace(mesh, 'CG', 2)

    # Define boundary conditions
    # Inlet velocity
    bcu_inlet = DirichletBC(V, inflow_profile(mesh, 2), surfaces, inlet_tag)
    # No slip
    bcu_wall1 = DirichletBC(V, Constant((0., 0.)), surfaces, wall_tag1)
    bcu_wall2 = DirichletBC(V, Constant((0., 0.)), surfaces, wall_tag2)
    bcu_wall3 = DirichletBC(V, Constant((0., 0.)), surfaces, wall_tag3)

    # Fixing outflow pressure
    bcp_outflow = DirichletBC(Q, Constant(0.), surfaces, outlet_tag)

    bcu_jet = DirichletBC(V, omega*g, surfaces, control_tag)

    bcu = [bcu_inlet, bcu_wall1, bcu_wall2, bcu_wall3, bcu_jet]

    bcp = [bcp_outflow]

    # Define trial and test functions
    u, v = TrialFunction(V), TestFunction(V)
    p, q = TrialFunction(Q), TestFunction(Q)

    # Define functions for solutions at previous and current time steps
    u_n = Function(V, annotate=annotate)
    u_ = Function(V, annotate=annotate)
    p_n = Function(Q, annotate=annotate)
    p_ = Function(Q, annotate=annotate)
    u_aux = Function(V)

    # Define expressions used in variational forms
    U = 0.5*(u_n + u)
    n = FacetNormal(mesh)
    f = Constant((0, 0))
    ix = Constant((1, 0))
    k = Constant(dtn)
    mu = Constant(mu)

    # Define variational problem for step 1

    F1 = rho*dot((u - u_n) / k, v)*dx \
        + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
        + inner(sigma(U, p_n, mu), epsilon(v))*dx \
        + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
        - dot(f, v)*dx
    a1 = lhs(F1)
    L1 = rhs(F1)
    # Define variational problem for step 2
    a2 = dot(nabla_grad(p), nabla_grad(q))*dx
    L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx
    # Define variational problem for step 3
    a3 = dot(u, v)*dx
    L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
    # Assemble matrices
    A1 = assemble(a1)
    A2 = assemble(a2)
    A3 = assemble(a3)
    # Apply boundary conditions to matrices
    [bc.apply(A1) for bc in bcu]
    [bc.apply(A2) for bc in bcp]
    [bc.apply(A3) for bc in bcu]

    As = [A1, A2, A3]

    # Create solvers
    solvers = list(map(lambda x: LUSolver(), range(3)))

    # Set matrices
    for s, A in zip(solvers, As):
        s.set_operator(A)

    u_start = project(u_init, V)
    p_start = project(p_init, Q)

    u_.assign(u_start, annotate=annotate)
    p_.assign(p_start, annotate=annotate)

    i = 1
    t = 0.0  # Initial time
    times = [t, ]

    import tqdm
    timer = tqdm.tqdm("Solving forward problem", total=num_steps)
    outfile_u = XDMFFile("output/u.xdmf")
    outfile_p = XDMFFile("output/p.xdmf")

    outfile_u.write(u_, t)
    outfile_p.write(p_, t)

    Jtemp = assemble(inner(curl(u_), curl(u_))*dx)
    Jlist = [Jtemp]

    while t <= T:

        timer.update(1)
        # Step 1: Tentative velocity step
        b1 = assemble(L1)
        [bc.apply(b1) for bc in bcu]
        solvers[0].solve(u_.vector(), b1, annotate=annotate)
        # Step 2: Pressure correction step
        b2 = assemble(L2)
        [bc.apply(b2) for bc in bcp]
        solvers[1].solve(p_.vector(), b2, annotate=annotate)
        # Step 3: Velocity correction step
        b3 = assemble(L3)
        solvers[2].solve(u_.vector(), b3, annotate=annotate)

        # Plot solution
        # print("Current value of u_ at time", t, "is:", u_.vector().get_local())

        # Update previous solution
        u_n.assign(u_, annotate=annotate)
        p_n.assign(p_, annotate=annotate)
        outfile_u.write(u_, t)
        outfile_p.write(p_, t)

        Jtemp = assemble(inner(curl(u_), curl(u_))*dx)
        Jlist.append(Jtemp)

        # # Update time
        t += float(dtn)

    outfile_u.close()
    outfile_p.close()

    return u_, Jlist, dtn


# Callback function for the optimizer
# Writes intermediate results to a logfile
def eval_cb(j, m):
    """ The callback function keeping a log """

    print("omega = %15.10e " % float(m))
    print("Functional = %15.10e " % j)


def optimize(dbg=False):
    """ The optimization routine """

    # Final time and time dicretization
    T = 0.5
    dtn = 0.1
    num_steps = int(T/dtn)

    # Define the control
    Omega = Constant(2e-1)
    source = Source(Omega, degree=2, name="source")

    # Execute first time to annotate the tape
    u, Jlist, dtn = forward(source, Omega, True)

    # Define the control
    control = Control(Omega)

    J = 0
    for i in range(1, len(Jlist)):
        J += 0.5*(Jlist[i-1] + Jlist[i])*float(dtn)

    # compute the gradient
    dJd0 = compute_gradient(J, control)
    print("gradient = ", float(dJd0))

    h = Constant(0.1)

    # Prepare the reduced functional
    reduced_functional = ReducedFunctional(J, control, eval_cb_post=eval_cb)

    # Taylor test for convergence rate
    conv_rate = taylor_test(reduced_functional, Omega, h)

    # Run the optimisation
    omega_opt = minimize(reduced_functional, method="L-BFGS-B",
                         tol=1.0e-12, options={"disp": True, "gtol": 1.0e-12})

    # Print the obtained optimal value for the controls
    print("omega = %f" % float(omega_opt))

    return omega_opt


optimize()
5

Dear Dokken,
I’ve tried to implement the dirichlet bc as you’ve suggested, but treating omega as a multiplicative constant seems to result in the Functional being independent of the control variable(omega).