Import of firedrake_adjoint leads to non-convergence

[hbuesing]

Dear all,

When I import firedrake_adjoint my Firedrake program does not convergence any more. I created the example below, which converges without firedrake_adjoint and does not with it. The issue has something to do with imposing bc0. Thank you for your help!

Henrik

# Code example

# coding: utf-8

'''
Minimum failing example 2p flow with dolfin-adjoint.
'''
from firedrake import *
from firedrake_adjoint import *

Nx=10
Ny=1
Nz=10

# Length in x- and y-direction. 
Lx=1.0
Ly=1.0
# Layer height of extruded mesh
Delta_z = 1.0/Nz
# Length in z-direction
Lz = Delta_z*Nz

# Construction of FV mesh
meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True)
mesh = ExtrudedMesh(meshbase, Nz, Delta_z)

def forward():
    horiz_elt = FiniteElement("DQ", quadrilateral, 0)
    vert_elt = FiniteElement("DG", interval, 0)
    elt = TensorProductElement(horiz_elt, vert_elt)
    DG1 = FunctionSpace(mesh, elt)
    DG2 = FunctionSpace(mesh, elt)

# Mixed Space
    W = DG1 * DG2

# Nonlinear problem --> no Trial functions
    u = Function(W)
    u_old = Function(W)
    pw, Sn = split(u)
    pw_old, Sn_old = split(u_old)

    v, q = TestFunctions(W)

# Define coordinate functions
#
    x_func_expr = Expression("x[0]")
    y_func_expr = Expression("x[1]")
    z_func_expr = Expression("x[2]")

    x_func = interpolate(x_func_expr, DG1)
    y_func = interpolate(y_func_expr, DG1)
    z_func = interpolate(z_func_expr, DG1)

# Define 1 for (1-Sn)
    one = Constant(1.0)

    uw_old = one-Sn_old
    uw = one-Sn
    un_old = Sn_old
    un = Sn

    dt = Constant(100.0)

# no capillarity atm. 
    F1 = inner((uw - uw_old)/dt,v)*dx
    F2 = inner((un - un_old)/dt,q)*dx
    F = F1 + F2
# Hydrostatic water pressure (for imposing pw Dirichlet BC)
    fun = Function(DG1)
    fun = (Lz -Delta_z/2 - z_func)*1e4+1e5

# Boundary Conditions
# pw BC (right)
    bc0 = DirichletBC(W.sub(0), fun, 2,method="geometric")
    t= Constant(dt)
    tend = Constant(1000.0)
    dtmax = 5000.0
    dtmin = 1e-3
    l = 1

    problem = NonlinearVariationalProblem(F,u,bcs=[bc0])
# Control PETSc options from the command line
    solver = NonlinearVariationalSolver(problem,options_prefix="",solver_paramaters={'ksp_type': 'preonly', 'pc_type': 'lu','snes_atol': '1e-6', 'snes_rtol': '1e-6','snes_stol': '1e-6','snes_type': 'newtonls', 'snes_linesearch_type': 'basic', 'snes_max_it': '20'})
    count = 0

    while (float(t) <= float(tend)):
       if COMM_WORLD.rank == 0:
         print("Timestep: ", l, "Cum. time: ", float(t), "Step: ", float(dt))
       try:
         solver.solve()
       except: 
         dt.assign(dt/2.0)
         u.assign(u_old)
         l = l+1
         if float(dt)<=dtmin:
           if COMM_WORLD.rank == 0 :
             print("Minimum step size reached. Aborting.")
           raise SystemExit(0)
       else:
         l = l+1
         if float(dt)<dtmax :
           dt.assign(dt*2)
         if float(dt)>=dtmax :
           dt.assign(dtmax)
         t.assign(t + dt)
         if ((float(t)>float(tend)) and count == 0) :
           t.assign(t - dt)
           dt.assign(tend -t)
           t.assign(t + dt)
           count += 1
         pw, Sn = split(u)
         u_old.assign(u)

forward()

[dham]

Hi @hbuesing, is this really minimal? E.g.:

• Do we really need lots of timesteps, or does this diverge after just one or two?
• Do we need all the terms, or does it still exhibit the problem with some removed?
• Do we need the variable timestep code, or does it still diverge with a small, fixed timestep?

[hbuesing]

I further reduced the example.

Thank you! Henrik

# coding: utf-8

'''
Minimum failing example 2p flow with dolfin-adjoint.
'''
from firedrake import *
from firedrake_adjoint import *

Nx=10
Ny=1
Nz=10

# Length in x- and y-direction. 
Lx=1.0
Ly=1.0
# Layer height of extruded mesh
Delta_z = 1.0/Nz
# Length in z-direction
Lz = Delta_z*Nz

# Construction of FV mesh
meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True)
mesh = ExtrudedMesh(meshbase, Nz, Delta_z)

def forward():
    horiz_elt = FiniteElement("DQ", quadrilateral, 0)
    vert_elt = FiniteElement("DG", interval, 0)
    elt = TensorProductElement(horiz_elt, vert_elt)
    DG = FunctionSpace(mesh, elt)

# Mixed Space
    W = DG * DG

# Nonlinear problem --> no Trial functions
    u = Function(W)
    u_old = Function(W)
    pw, Sn = split(u)
    pw_old, Sn_old = split(u_old)

    v, q = TestFunctions(W)

# Define coordinate functions
#
    z_func_expr = Expression("x[2]")
    z_func = interpolate(z_func_expr, DG)

# no capillarity atm. 
    F = inner(Sn - Sn_old,v)*dx

# Hydrostatic water pressure (for imposing pw Dirichlet BC)
    fun = Function(DG)
    fun = z_func*1e1

# Boundary Conditions
# pw BC (right)
    bc0 = DirichletBC(W.sub(0), fun, 2,method="geometric")

    problem = NonlinearVariationalProblem(F,u,bcs=[bc0])
# Control PETSc options from the command line
    solver = NonlinearVariationalSolver(problem,options_prefix="",solver_paramaters={'ksp_type': 'preonly', 'pc_type': 'lu','snes_atol': '1e-6', 'snes_rtol': '1e-6','snes_stol': '1e-6','snes_type': 'newtonls', 'snes_linesearch_type': 'basic', 'snes_max_it': '20'})
    count = 0
    try:
      solver.solve()
    except:
      print("No convergence!")
    else:
      print("Convergence!")

forward()

[dham]

OK. I agree that that's pretty minimal. I will try to find time imminently to walk through this. I wonder if it's the geometric BCs. That's a rarely used feature which we might have somehow dropped in the adjoint code...

[dham]

I can reproduce this. Now trying to find out what's actually happening...

[dham]

OK, one convenient flight from Munich to London later, I can see what's happening. The try clause is obscuring the actual error. Take that off and you find that dolfin-adjoint is raising an exception because it doesn't know how to annotate the assignment of UFL expressions. This is happening when the Dirichlet boundary conditions are being applied. The problematic line is therefore actually:

fun = z_func*1e1

This doesn't actually set the values of the Function fun to zfunc*1e1, instead it sets the variable fun to the expression zfunc*1e1. To work around this, instead use:

fun.interpolate(z_func*1e1)

As an aside, the code provided also misspells solver_parameters so the solver options are not actually applied.

[hbuesing]

Dear David,

Thank you very much for sorting this out! Best regards and "Happy Easter"!

Henrik

-- Dipl.-Math. Henrik Büsing Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University

Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889

http://www.eonerc.rwth-aachen.de/GGE hbuesing@eonerc.rwth-aachen.de

---

Von: David A. Ham [notifications@github.com] Gesendet: Donnerstag, 13. April 2017 21:24 An: firedrakeproject/firedrake Cc: Buesing, Henrik; Mention Betreff: Re: [firedrakeproject/firedrake] Import of firedrake_adjoint leads to non-convergence (#1042)

OK, one convenient flight from Munich to London later, I can see what's happening. The try clause is obscuring the actual error. Take that off and you find that dolfin-adjoint is raising an exception because it doesn't know how to annotate the assignment of UFL expressions. This is happening when the Dirichlet boundary conditions are being applied. The problematic line is therefore actually:

fun = z_func*1e1

This doesn't actually set the values of the Function fun to zfunc1e1, instead it sets the variable fun to the expression zfunc1e1. To work around this, instead use:

fun.interpolate(z_func*1e1)

As an aside, the code provided also misspells solver_parameters so the solver options are not actually applied.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/firedrakeproject/firedrake/issues/1042#issuecomment-293998062, or mute the threadhttps://github.com/notifications/unsubscribe-auth/APAHzryEvPZLmQKgzEVoYDF-UhdZtlqeks5rvnaJgaJpZM4M43Cv.